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Coordinate  Conversion  Technique 
for  OTH  Backscatter  Radar 


I.  INTRODUCTION 

The  WIMP  Program  is  a complex  algorithm  which  purports  to  permit  accurate 
conversion  from  radar  target  range  coordinates  to  geographic  target  coordinates  in 
OTH  Backscatter  ITF  radar  applications;  the  complexity  is  fundamentally  due  to  a 
deeply  developed  iterative  and  reiterative  algorithm  that  ultimately  involves  the 
user.  The  program  consists  of  three  main  components-. 

(a)  A set  of  subprograms  which  construct  a three-dimensional  three-layer 
ionospheric  model. 

(b)  A set  of  subprograms  which  simulate  IIF  radio  wave  propagation  in  the 
ionosphere  constructed  in  (a). 

(c)  One  of  three  driving  programs  controlling  the  calling  sequence  to  the  radio 
wave  propagation  simulator  in  (b > ; these  driving  programs  perform  the  following 
tasks:  (1)  From  a given  ionosphere,  generate  the  leading  edge  of  a simulated  ob- 
lique ionogram  including  scale  factors  to  the  Fg-layer  parameters  to  force  agree- 
ment with  a given  vertical  ionogram;  (2)  from  a comparison  of  the  simulation  in 

(1)  to  a given  oblique  ionogram  generate  range  gradient  factors  to  apply  to  ^Fg  and 
MfSOOOlFg  to  force  agreement;  (3)  from  the  final  ionosphere  generated  in  (1)  and 

(2) ,  simulate  predicted  radio  frequency  propagation  paths,  from  which  range  vs 
group  path  functions  may  be  tabulated. 

(Received  for  publication  24  May  1977' 
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We  first  describe  in  some  detail  the  WIMP  three-dimensional  three-layer 
ionospheric  model,  followed  by  a general  description  of  the  radio  wave  propagation 
simulation  subprogram  (ray-tracing  algorithm).  A functional  description  is  gi\  m 
of  the  three  driving  programs;  the  operating  instructions  for  their  use  are  in- 
cluded next. 

Program  block  diagrams,  flow  charts,  and  program  listings  are  included  in 
Appendixes  A,  B,  and  C. 

Comparisons  have  been  made  of  the  results  of  the  WIMP  ray  tracings  with 
techniques  developed  at  It  A DC . The  description  and  results  or  these  studies  are 
reported  in  Sections  19  through  21. 


2.  THU  WIMP  Til  IONOSPHERIC  MODEL 

A complete  three-dimensional  ionospheric  model  generally  consists  of  three 
fundamental  aspects;  (1)  How  many  distinct  layers  are  to  be  included,  (2)  How 
shall  the  layer  parameters  be  modelled,  (3)  How  shall  the  electron  density  profile 
be  erected  from  the  modelled  parameters?  A collateral  question,  equally  important 
is,  what  ray-tracing  technique  is  to  be  used,  that  is,  which  aspect  dominates  the 
simulation  of  I1F  propagation?  To  put  the  following  discussion  in  proper  perspective 
the  latter  question  will  be  considered  here  in  general,  reserving  a more  detailed 
treatment  to  a latter  section. 

The  WIMP  ray-tracing  algorithm,  SUBROUTINE  TRISL,  is  a control  point 

technique  based  on  layer  parameters  as  in  the  ITS-7  8 program*  and  the  DEVAN 
2 3 

program.  ’ The  treatment  of  gradients  of  electron  density,  however,  differs; 
ITS-78  does  not  consider  gradients  explicitly,  DEVAN  constructs  an  equivalent 
tilted  reflecting  surface  from  variations  in  virtual  height,  while  WIMP  constructs 
the  tilt  of  the  reflecting  layer  from  the  gradients  of  electron  density.  The  predic- 
tions of  the  WIMP  model  are  primarily  sensitive  to  the  layer  parameter  predictions, 
and  only  weakly  sensitive  to  the  electron  density  profile  model.  This  latter  de- 
pendence only  enters  the  calculation  at  the  reflection  point  (which  generally  is  within 
the  F 2 layer),  and  involves  the  horizontal  gradients  of  the  electron  density  which 
are  determined  primarily  by  the  gradients  of  the  Fg -layer  parameters. 

The  WIMP  model  ray  tracing  considers  three  distinct  layers,  E,  Fj,  and  Fgi 
an  optional  D-layer  tail  may  be  appended  to  the  E layer.  The  E-layer  critical 

1.  Barghausen,  A.  F.  , et  al  (1969)  Predicting  Long-Term  Operational  Parameters 

of  High-Frequency  Sky-Wave  Telecommunication  Systems,  ESbA  'Technical 
Report  FRL  llO-lTS-78,  Institute  for  Telecommunication  Sciences, 

Boulder,  CO. 

2.  Beckwith,  R.I.,  Bailey,  A.  D.  , and  Rao,  N.  N.  (1972)  An  Investigation  of  Direc- 

lional  Propagation  Effects  in  High-Frequency  Radio  Source  Location,  RRL 
Publication  No.  409. 

3.  Beckwith,  R.  I.  (1973)  A Computer  Program  for  the  Rapid  Prediction  of  Angles 

of  Arrival  of  IIF  Radio  Waves,  RRL  Publication  No.  442. 
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frequency  is  determined  from  the  local  solar  zenith  angle.  The  model  employed  for 
the  latter  is  sufficiently  unique  to  describe  in  some  detail  below.  The  Fj  critical 
frequency  is  25  percent  larger  than  f E with  an  additional  0.5  MHz. 

The  ITS  model  provides  predictions  for  the  F ^ critical  frequency  and  M13000) 
factor.  These  parameters  may  be  modified  by  ad  hoc  scale  factors  to  foi  cc  agree- 
ment with  a given  vertical  ionogram,  and  ad  hoc  gradient  factors  to  force  agree- 
ment with  a given  oblique  ionogram.  The  Eg  semithickness  is  determined  from  a 
linear  sunspot  model.  All  parameters  thus  far  mentioned  are  used  to  generate  the 
F2  -layer  height. 

The  Fj  layer  is  fitted  between  the  E and  layers  such  that  the  F^  bottom 
coincides  with  the  E-layer  height  and  the  F^ -layer  height  coincides  with  the  F2- 
layer  bottom. 

A more  detailed  description  of  the  layer  parameters  and  electron  density  model 
follows. 


3.  MODEL  FOR  SOLAR  ZENITH  ANGLE 

The  model  for  f E requires  the  solar  zenith  angle.  From  the  law  of  cosines  of 
a spherical  triangle,  this  quantity  is  given  by: 

cos  Z = sin  6 sin  0 + cos  5 cos  9 cos  (<f>  - 0 ) , (1) 

s 

where 


Z is  the  solar  zenith  angle, 

6 is  the  solar  declination  (=latitude  of  the  subsolar  point), 
bs  is  the  east  longitude  of  the  subpolar  point, 

6 is  the  geographic  latitude  of  the  field  point, 
d>  is  the  east  longitude  of  the  subsolar  point. 

The  field  point  (0,  <$)  is  given,  and  the  subsolar  point  (6,  4>  ) is  required. 

The  subsolar  point  may  be  determined  from  the  solution  of  the  Kepler  problem 
using  mean  elements  from  the  American  Fphemeris  and  Naufical  Almanac.  Such  a 
solution,  to  second  order  in  the  eccentricity  of  the  earth's  orbit  a<  d neglecting  lunar 
perturbations,  produces  an  ephemeris  accurate  to  within  a nautical  mile. 

Alternatively,  a simplified  model  for  the  solar  ephemeris  may  be  constructed, 
accurate  to  one  degree  of  arc,  as  follows: 


6 = E sin  A,  A = 27r(d-80)/365.  25  , 
0S  = 7T  - UT  , 


(2) 
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where,  F is  t he  obliquity  of  the  earth's  equator  (=23.°45),  A is  the  mean  longitude 
of  the  sun  measured  in  the  ecliptic  counterclockwise  from  the  first  point  of  Aries, 

FT  is  the  universal  time,  and  d the  day  of  year. 

Clearly,  models  between  these  two  may  be  constructed;  in  particular,  the 
annual  variation  of  the  equation  of  time  may  be  modelled  as  a correction  to  the 
subsolar  longitude.  (The  equation  of  time  is  the  correction  to  be  added  to  "sundial 
time"  to  reduce  it  to  local  mean  time.  ) Such  an  approach  is  used  in  the  WIMP  pro- 
gram, with  the  following  algorithm  implemented  in  subroutine  PTHP  each  time  it  is 
called.  The  solar  declination  ("zenith  correction  angle  for  day  of  year")  is  approxi- 
mated bv: 


6'  = E sin  (2jr(d-82.  5)/364)  , (3) 

The  equation  of  time  in  hours  is  modelled  by: 

X.  = aYv  [A  1/(B1  + X j)]  - [A2/(B2  + X2»  . (4) 

where 

a 0. 0235  hours; 

Y = * 1.  , X2  = d + 21.  , 0.0  < d < 162.5  ; 

Y -1.  , X2  = 344. 5 - d,  162. 5 * d < 344. 5 ; 

Y = + 1.  , X2  - d - 344.  5,  344.  5 £ d < 366.  5 ; 

Xj  = (X2/ 45.  )3 ; 

A j =64.  , A2  = 240.  , Bj  = 4.  , B2  = 15. 


The  following  angle  is  defined,  being  converted  to  radians: 


D - 0'  + X.  - 0 , 
J 


(5) 


where 


0’  = 71  - UT,  UT  > v , 


0'  = 7T  t UT,  IJT  < 77 


The  following  approximation  to  the  solar  zenith  angle  is  constructed: 
2 


Z 


(D*cos(0.  6(0  + 6')) 


2 * (0  - 6')2  , 


(6) 
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and  the  cosine  of  the  solar  zenith  angle  is  approximated  by: 


cos  Z = cos  (0.932  Z&)  . (7) 

In  the  above  description,  angles  are  considered  in  radians,  including  UT; 
furthermore,  longitudes  are  conventionally  taken  to  be  positive  east.  In  the  coding 
of  PTHP,  angles  are  conventionally  in  degrees,  UT  in  hours,  and  longitudes  are 
positive  west. 

t.  MODEL  FOR  THE  E-LAYER  PARAMETERS 

The  E-layer  critical  frequency  fg,  layer  height  h^,,  semithickness  y^,,  and 
bottom  altitude  h^  are  required.  If  the  simple  parabolic  E -layer  model  is  selected, 

h,  = 100  km, 

b 

h^  = 115  km,  (8a) 

h, 

yE  = 15  km. 

If,  on  the  other  hand,  the  "parabolic  E-  with  D-layer  tail"  electron  density  model 
is  selected, 

hb  = 100  -4.FlF2, 

hE  = 115  km,  (8b) 

yE  hE  ‘ \ ' 
where 


Fj  = 1 - (1  -R/3000F 

R < 6000  km 

o 

II 

R F 6000  km 

F9  = cos(A) 

A < n/2  : 

CO 

11 

o 

A > tt/2  ; 

R = DIST  (6,  0,  76.  0°  N, 

102.0°  W); 

A = (i>+^-|t-7r|-^7i; 

and  (6,  0)  is  the  geographic  north  latitude  and  east  longitude  of  the  field  point; 
DIST  is  the  great  circle  distance  (in  km)  between  the  two  points  specified  in  the 
argument  list. 
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The  critical  frequency  fj,  is  modelled  from  the  sunspot  number  and  the 
cosine  of  the  solar  zenith  angle  cos  Z [as  modelled  by  Kq.  (7)].  bet 

1 

C . 1(1  i cos  7.)2,  S,  = (1  t 0.004  N )6 7  . 

14  1 s 

f 0.75S.C.,  f 3.  17  S.  (cos  Z)a  , (9a) 

n 1 T x 1 

1 

r . - (f  2 t f 2)7,  a = 1/3.775  . 
d x n 

Then,  at  night  (cos  Z negative)  f\p  is  given  by  f ; during  the  day  (cos  Z positive)  fr<, 

is  given  by  f^,  unless  exceeds  f (I^Fg  as  modelled  by  the  ITS-78  model)  in  which 

case  f is  taken.  Furthermore,  f„  is  constrained  not  to  exceed 
n E 

r_  = (0.  99  f_  - 0.  5)/ 1.  26  . (9b) 

E max  2 

5.  MODEL  FOR  THE  F, -LAYER  CRITICAL  FREQUENCY 

The  Fj-layer  critical  frequency  f^  is  given  by 
fj  = 0.5  + 1.  26  rE  , 

and  is  constrained  [ Eq.  (9b)]  not  to  exceed  0.  99  f . 

6.  MODEL  FOR  f()F2  AND  M(3000)F2 

These  two  parameters,  and  M^,  are  given  by  the  ITS-78  parameter  predic- 
tion model.  This  model  is  sufficiently  well  known  that  it  will  not  be  descr  ibed  in 
detail  herein.  Generally,  the  model  consists  in  reducing  ...  set  of  coefficients  for 
each  modelled  parameter  to  a given  phase  of  the  solar  cycle,  season  of  the  year, 
universal  time,  and  geographic  position.  The  solar  cycle  is  represented  by  the 
sunspot  number,  and  the  parameters  are  fitted  to  a first  or  second  order  poly- 
nomial. The  seasonal  variation  of  the  various  parameters  is  represented  by  differ- 
ent sets  of  coefficients  for  each  month,  except  for  ^F^.  The  latter  is  represented 
as  a Fourier  sum  (9  terms)  in  the  day  of  the  year.  The  diurnal  variation  is  repre- 
sented by  a Fourier  expansion  (up  to  13  terms)  and  the  global  variation  is  represented 
by  a (modified)  spherical  harmonic  expansion.  Parameters  modelled  are 
M(3000)F9,  h.  F„,  f E,  upper  decile  f E , median  f E , and  the  lower  decile 
f E . Only  the  first  two  parameters  are  used  herein. 


(10) 

(10) 
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7.  MODKI.  FOR  TIIK  F2-1.AYKR  SKMITIIICKNF.SS 

The  F.. -layer  semithiekness  is  given  by 

v„  = F (56  1 0.489  N ) . (11) 

2 y s 

The  model  allows  lor  a correction  scale  factor  F , which  is  initialized  to  unity;  no 

y J 

provision  is  made  for  its  redefinition. 


MODKI.  FOR  TIIK  K2-!.AYKR  IIKIGIIT 

The  model  employed  for  the  construction  of  the  Fg -layer  height  uses  all  of  the 
parameters  thus  far  described;  f^.,  fj,  f^,  M ^ . Yg.  h^.  The  following  fre- 
quencies are  defined; 

r,  - 

r f2  - 0.  01  f2  y2, 

f = max  (f1  , f"), 

q q q 

f3  r2M3* 

f0  ' 4 fq  - 3 f2’ 

f = f + 3.78  f„  + 1.5, 


1.  26  fE  1 0.  5, 


(12) 


o 


8.  82  f + 3.  5, 


r 


1.  26  f„  - 0.  5 = f = f,, 
h o 1 


o 

8. 00  fE, 

ro  t 7.78  fE  1 1.5, 

f - 0.  22  f„  r 1.5, 
o L 

f9  f fn  * 

2 q 

r9  - r . 

2 q 


The  following  parameters  are  defined  (these  are  stored  for  later  use  in  the 
construction  of  the  IVHSOOOlFg  correction  factor; 


Il3  = 70  t 3100  (f„/f9)  , 


q 3 


F = f / f 

q q 2- 


(13) 
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K,  <ra/rb)  ioH  <fc/r„>. 

Z1  ’ yE  VV  '°e  W’ 
z2  iv2'Vr2"oe'ri,fj'- 

'I'he  F'^ -layer  height  is  given  by: 

h2  8,0  <H3  " Z\  " Z2  _ 'V/P3  + y2  + hE  ‘ 

9.  MODEL  FOK  TIIK  F,-LAYER  HEIGHT  AND  SEMITHICKNESS 

The  Fj-layer  height  hj  and  semithickness  y^  are  given  by: 

hl  h2  - y2’ 


(14) 


(15) 


10.  MODEL  FOK  F*>- LAYER  CORRECTION  FACTORS  - 
VERTICAL  10  NOG  RAM 


The  following  information  had  been  extracted  from  a vertical  incidence  ionogram 
as  input  to  the  program  from  which  correction  factors  to  the  Fg  layer  parameters 
are  to  be  derived: 

t — universal  lime  of  ionogram, 

L — geographic  north  latitude  of  sounding  station, 

W — geographic  west  longitude  of  sounding  station. 


Fj  — f F'j  from  ionogram 
P2 


— fQ lv>2  from  ionogram. 


F 


F’ 

q 


IT 


D 


1 


4 (F2 
1 


Fl>’ 


" F2-T(F2-F1)' 

— Fg  virtual  height  at  frequency  F , 

— Fg  virtual  height  at  frequency  F^, 

— sunspot  number, 

— day  of  year. 


The  ionospheric  parameters  at  the  given  event  are  constructed  from  the  param- 
eter prediction  model;  let  these  parameters  be  represented  by  the  following  notation: 


14 


hj..,  hj,  h 

■V  -Yr  y 


2 

2 


\ 


— critical  frequencies, 

— layer  heights, 

— semitheikness, 

— bottom  of  ionosphere, 


M, 


M(3000)F2, 


f — "quarter  point  frequency",  Eq.  (14), 

Zj,  Z^,  Pg  — Eq.  (15)  parameters. 

Define  the  following  frequencies: 


f' 

q 


f 

a 


f 


cl 


= fn 


ff  f2 


f = f 
h 2 


f. 

l 


+ 3.78  f + 
+ 8.  82  f + 


1.5, 
3.  5, 


7.7  8 f + 1.5 

- 0.  22  f + 1.  5, 

+ P . 

q 

- p . 

q 


Define  the  following  quantities: 


Zj  = (yE/8)*(fa/fE)*  log  (ff/fh), 

Z2  ^2,2)*%,t2)*  ,0*  (fj/fj). 

P3  = (f. a/rj)*  log  (fc/fd), 

H3  <P3/P3)*(Hv-Z'l-Z21-\}+Zl  + Z2,hb’ 

M'  (F/fJ*( 3100/H’  -70))?. 

j C]  z J 

Then  the  correction  factors  Rj.  and  purportitig  to  correct  predicted  ^Fg 
and  M(3000)F2  parameters  are  given  by: 


Rf  I*  2^2' 
Rm  ^/M3  • 


(16) 
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1 1 . M0DK1,  FOR  Fv-LAYKR  CORRKCTION  FACTORS  - 
OHLIQUK  IONOGRAM 


Factors  IrR  and  MR  are  calculated  in  one  run  of  the  OBLFACT  program,  and 
are  used  as  input  for  die  next  run  of  OBLFACT  to  be  stored  in  elements  3 and  5 
of/PERTC/.  The  reference  point  is  also  input  directly  into  elements  1 and  2 of 
/PERTC/;  the  remaining  elements  are  established  by  DATA  statements.  In  the 
following  description,  the  elements  of  /PERTC/  are  referenced  by  the  elements  of 
(Cji  i=  1,  16).  The  field  point  north  latitude  and  west  longitude  (N,  W)  are  arguments 
of  PERTC,  the  algorithm  for  the  calculation  of  the  oblique  correction  factors.  The 
reference  latitude  and  longitude,  C^  and  C^,  are  alternatively  referenced  as  and 
W , respectively. 

The  algorithm  commences  by  computing  the  distance  and  azimuth  to  the  field 
point  (N,  W)  from  the  reference  point  (N^,  WM;  let  these  two  quantities  be  D and  Z , 
respectively. 

Define  the  following  quantities: 


Zj  (ir/180)*(C13  - C12  - 90), 

S esc  (Z  j), 

Y 3600  Cg  C j j S, 

W = c1Qs 

Z2  (W180)*((D  t Y)/360  W), 

X sin  (Zg). 

The  following  correction  factors  are  then  calculated: 


R =(DC.X)*(DC.D/C.)*(DC  Z)  (17) 

a l j k m 

where  the  f F^  and  M(3000)F2  corrections  are  specified  by  the  index  set  (a;  i,  j,  k,  m) 
being  assigned  values  (f;7,  3,  4.  14)  and  (m;f),  5 6,  15),  respectively.  These  factors 
are  then  applied  to  the  appropriate  1?2 -level  parameter. 


12.  Kl,  KM  HINTS  OF  COMMON  HLOCK/l’KRTC/ 

This  common  block  is  linked  to  the  INPUT  file  via  NAMELIST  inputs  in  the  main 
programs  MAIN  and  FTDBLB;  OBLFACT  links  elements  of  NAMELIST  input  (FRO 
and  MRO)  as  arguments  of  subroutines  SKIP  and/or  PEAK;  the  latter  two  routines 
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then  define  elements  3 and  5 of  /PERTC / equal  to  FRO  and  MRO,  respectively. 
OBI, FACT  also  links  the  first  two  elements  NFREF  and  WLREF  to  the  input 
NAMELIST.  All  elements  are  assigned  default  values  in  a BLOCK  DATA  program; 
presumably,  these  defaults  may  be  assigned  different  values  by  replacing  the 
BLOCK  DATA  program.  The  elements  of  / PE RTC/  are  used  by  SUBROUTINE 
PERTC  which  computes  additional  corrections  to  f F„  and  M(3000)F„  based  on 

O fa  Cj 

information  extracted  from  an  oblique  ionogram.  The  reference  point  (NLREF, 
WLREF)  is  the  same  reference  point  used  in  the  construction  of  the  updating  cor- 
rection factors  to  f F„  and  M(3000)Fo  effected  in  SUBROUTINE  FACTO. 

O u C. 


Table  1.  The  Elements  of  /PERTC/  with  the  Mnenomic  Names  Used  in  PERTC 
and  the  Default  Values  Assigned  in  the  BLOCK  DATA  Program 


Element 

Mnenomic 

Default 

Element 

Nmenomic 

Default 

1 

NLREF 

2 

WLREF 

3 

FR 

0.  0 

4 

FD 

1000.  0 

5 

MR 

0.  0 

6 

MD 

1000.  0 

7 

FZ 

0.  0 

8 

MZ 

0.  0 

9 

DUT 

0. 0792 

10 

DW 

14  8.  5 

11 

X 

0.  225 

12 

GAMMA 

180.  0 

13 

AZM 

38. 808 

14 

DFAC 

0.  0 

15 

DMAC 

0.  0 

16 

AZREF 

0.  0 

13.  ELECTRON  DENSITY  MODEL 


A continuous  electron  density  profile  is  constructed  from  the  predicted  ionos- 
pheric parameters.  Slope  discontinuities  may  occur  in  the  model  at  the  critical 
points  of  the  profile.  The  model  allows  for  either  a parabolic  model  for  the  lower 
E-layer  (below  the  critical  altitude)  or  an  E-layer  with  a D-layer  tail;  these  two 
models  will  be  described  first. 

In  the  parabolic  E-layer  model,  the  base  altitude,  layer  height,  and  semi- 
thickness are  modelled  to  be  ICO,  115,  and  15  km,  respectively.  Below  the  E-laycr 
height  hp,  the  electron  density  N^  is  given  by: 


N = 12400  f 2 (1  - Y2),  h,  < h < h„; 
e E b E 

N =0,  h < h,  , 

e b 


(18) 
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where 


Y - (hp  - h)/yj,;  h is  the  altitude;  hp,  y j , , and  fp  are  the  E-laye r height, 
semithickness,  and  critical  frequency. 

In  the  E-layer/D-layer  tail  model  the  semithickness  may  increase  by  up  to 
4 km,  hp  will  lower  correspondingly,  and  lip  remains  at  115  km.  First,  an  electron 
density  Neo  is  calculated  by  Eq.  (20)  with  Y,  however,  constructed  as  follows; 

V’  = 2 (hE  - h)/yp, , , 

x = 1.  5 | 1.7/(1. 7 + Y')]  ^ , 

z = [ 1 + 0.  1 (Y‘  - i)  | Y ' - l|  ] /[  1 - 0.  1/(Y'  + ll2]  , 

Y"=  2 - log  (1  + 1.72  (Y')x), 

Y = 1,  Y"  negative, 

Y = 1 - (Y")z,  Y"  not  negative. 

Then 

N = 12400  f 2 (1  - Y2)  . (10) 

eo  E 

Calculate  the  correction  for  the  D-layer  tail  as  follows: 

C = cos  (0.  85  X),  where  X is  the  solar  zenith  angle; 

= 0,  C negative;  j (20) 

Nd  = (1  + 0.004  ISM  CTHl  + 81  W4); 

where 

W = (65  - h)/y„,  and  N is  the  sunspot  number. 

E s 

The  electron  density  is  finally  given  by: 

N = N + N , . (21) 

e eo  d 

Between  hp  and  II^,  the  plasma  frequency  F is  taken  as  the  positive  solution 
of  a quadratic  equation: 

aF  2 + bF  + c = 0 (22) 

P P 
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and  the  electron  density  is  token  as 


N = 12400  (F  )2  . 
e p 


The  coefficients  a,  b,  and  c are  calculated  as  follows: 
First,  define  three  reference  heights  h^,  h^,  hc: 


1 

o 7 

V’,  'h  " - VV  > ’ 

hb  = h2  " y2  ^ ’ 


(23) 


For  h between  hr,  and  h , define  coefficients  B,,  B_,  Bn,  B.  as  follows: 

c 1 z o 4 


B1  2 (ha  ' hE  + ^3  ' hi)/,(fi  ‘ fE>' 

B2  = (hi'  " ■?  (fl(ha  " hE)  + fE(hb  " hl^(fl  " fE'’ 
B3  = y 1 /f  1 ’ 

B4  = yi  * 


For  h between  hc  and  hg,  the  B. 's  are  given  by: 

B1  =I,hb  ' hl1,,t2  ' V’  B2"  VVV 

B3  = y2^f2  ' B4  = y2  ’ 


1 ie  coefficients  a,  b,  c of  the  quadratic  Eq.  (22)  are  defined  as: 
2 2 

a : Bj  + B3  , 
b = 2 Bj  (B2  - h), 
c = (B2  - h)2  - B42  . 


The  topside  model  (h  greater  than  h^)  is  a modified  Chapman  profile.  Let 
Y = 2 (h  - hn) /y„;  then, 

Lt  C. 


Ng  = 1 2400 1 (f 2> 2 exp  [1  - Y - exp  (-  Y)] 

+ (0.  5 + 0.  1 f2)  h2  (h  - h2)/h2}  . 


(24) 
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I t.  MODI' I.  I'Olt  Till-:  (iUAIHKNT  OF  TIIK  KLKCTRON  DENSITY 


Any  ray-tracing  technique  will  require  the  gradient  of  the  electron  density  at  a 
field  point  (r,  0 , 0).  Clearly,  the  complications  of  the  model  described  above  for 
the  electron  density  mitigate  against  any  attempt  to  derive  analytic  expressions  for 
the  gradient,  even  neglecting  variations  in  the  ionospheric  parameters.  Further- 
more, should  one  desire  to  replace  the  electron  density  model  with  some  other, 
either  more  or  less  complicated,  a general  numerical  algorithm  for  the  gradient 
is  desired. 

The  gradient  model  employed  in  the  W IMP  program  places  the  field  point  at 
the  center  of  a cubical  volume,  and  computes  the  electron  density  at  each  of  the  nine 
points  so  defined.  The  coordinates  of  these  nine  points  are  listed  in  the  following 
table. 


Table  2.  Coordinates  of  Points  Used  for  Gradient  Calculation 


Point 

Radial  Coordinate 

Theta 

Phi 

1 

r 

9 

0 

2 

r + D 

r 

e , D0 

0+D0 

3 

r + D 

r 

6 -Dd 

*+D0 

4 

r + D 

r 

0 + D0 

^-D0 

5 

r+  D 

r 

e-De 

0 ~ D0 

6 

r - D 

r 

9 + De 

0 + »0 

7 

r - D 

r 

e-De 

0 + D0 

8 

r - D 

e + d„ 

0 - Da 

r 

9 

0 

Q 

r - D 

r 

e-uQ 

«'D0 

(r,  0 , 0)  are  the  coordinates  of  the  field  point, 
(D^,  D^,  D ) are  the  grid  spacings. 


D , - D / r sin  0 . 

0 r 

Given  a function  f(x)  evaluated  at  three  points,  x , x = x + d , and  x = x - dx, 

o + o x — o 

then  the  best  estimate  of  the  derivative  af  f at  xq  is  given  numerically  by: 

P = -l  (f  - f ) /dx2  (251 

o ^ + 

where  f , f , and  f represent  f(x  ),  f(x  1,  and  f(x  1.  The  best  estimate  of  the 
— o + - o + 

second  derivative  is  given  by: 


20 


(26) 


f"  = 4 (f  - 2f  + f )/dx2  . 

Unfortunately,  the  defined  set  of  points  does  not  contain  quite  the  correct  subset  to 
implement  the  numerical  model  for  the  derivative  specified  by  Fq.  (25).  However, 
appropriate  partial  derivatives  may  be  evaluated  at  the  center  of  each  side  of  the 
defined  cube.  Each  coordinate  then  has  four  approximations  to  the  appropriate 
partial  derivative  of  the  electron  density;  for  instance,  the  r partial  derivative  will 
use  pairs  (2,  G),  (3,  7),  (4,  8),  and  (5,  9). 

The  algorithm  adopted  to  represent  the  partial  derivatives  is  the  average  of  the 
four  approximate  values. 

The  value  chosen  for  the  grid  space  is  1 / 3-1/2  km. 


15.  THE  WIMP  3-1)  K W-TRACIISG  TECHNIQUE 

The  ray-tracing  model  employed  in  the  WIMP  program  depends  upon  the  appli- 
cation of  Bouger's  rule  aid  parabolic  electron  density  models  in  multiple  layers; 
the  gradient  of  the  elect’ on  density  at  the  reflection  height  is  employed  to  determine 
the  tilt  of  the  reflecting  layer. 

The  program  alio  </  • for  refraction  and/or  reflecting  in  the  E,  Fj,  and  Fg 
layers.  The  F layer  is  inc  hed  between  the  E and  F0  such  that  the  bottom  of  the 
Fj  layer  is  at  the  E-layer  ueight  and  the  F^-layer  height  is  at  the  bottom  of  the  Fg 
layer. 

The  ray  tracing  is  accomplished  in  SUBROUTINE  TRISL. 

15.1  IJouger’s  Rule 

Bouger's  rule  states  that  in  a refracting  medium  with  spherical  symmetry 
electromagnetic  radiation  propagates  in  such  a way  that  the  quantity  jur  cos  /3  re- 
mains constant  where  /i  is  the  local  index  of  refraction,  r is  the  distance  to  the  field 
point  from  the  center  of  symmetry,  and  /3  is  the  local  elevation  angle. 

15.2  Refraction  Through  a Layer 

Refraction  through  a layer  is  depicted  in  Figure  1.  The  group  path  correction 
AP'  phase  path  correction  AP,  and  range  correction  A are  given  by 

Ap,=  (ihwr){iul°g  (£M-)- 11- 

AP  = |l--^sin2/3c(l+l/u2)j>AP'--iysin/3c/u2,  (27) 

A = AP'  cos  B /(R  + h ), 

' c o m 
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where  u = f sin  ji  I f and  f is  the  operating  frequency,  f the  critical  frequency, 
Hq  the  radius  of  the  earth,  hm  the  layer  height,  y the  layer  semithickness,  and  ji 
the  local  elevation  angle  of  the  unrefracted  ray  at  the  layer  height. 


Figure  1.  Refraction  Through  Layer 


The  group  path  length  P'  and  phase  path  length  P upon  exiting  the  layer  are 
given  by 

P'  = D + AP'  , 
o 

(28) 

P = D + AP, 
o 

where  D is  the  distance  P'  P"  in  Figure  1. 

15.3  Reflection  from  ;i  Layer 

Reflection  from  a layer  is  depicted  in  Figure  2.  The  group  path  corrections 
A P are  given  by 
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AP'  = 


s 1 n p ~ u - 


AP  = ■j  y sin  ft  r + l1  ~ \ sin2  /3  r (1+  1 /u2)  {ap1  , 


(29) 


where  u = f sin  (1  ^/f  , and  is  the  local  elevation  angle  at  the  virtual  height  P". 
These  corrections  apply  for  the  upward  trace  to  the  reflection  point;  equal  correc- 
tions obtain  for  the  downward  trace. 


15.4  Tilted  Reflection  Layer 

The  gradient  of  the  electron  density  at  the  ray  reflection  point  is  evaluated  and 
employed  to  construct  an  effective  tilted  reflecting  layer.  Underlying  refracting 
layers  are  considered  concentric  with  the  earth;  however,  gradients  in  these  layers 
are  taken  into  account  in  that  the  upward  and  downward  legs  of  the  ray  trace  through 
the  underlying  layers  are  spatially  separated. 

The  introduction  of  a tilted  reflecting  layer,  of  course,  generates  considerable 
complication  in  the  ray  trace  in  that  the  ray  direction,  layer  normal,  and  radius 
vector  are  no  longer  coplanar.  The  progr  am  accounts  for  these  bookkeeping  com- 
plications by  setting  up  orthonormal  triads  as  the  ray  attacks  successive  layers. 
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The  y axis  is  taken  parallel  to  the  ray  direction.  The  yz  plane  contains  the  center 
of  curvature  of  the  tilted  reflecting  layer  (or  the  center  of  the  earth  in  a refracting 
layer). 

15.5  Iteration  Scheme 

In  each  hop  of  the  ray  trace,  the  various  calculations  are  iterated  to  refine  the 
location  of  the  profiles  at  which  refraction  through  a layer  and  reflection  from  a 
layer  occur.  In  the  reflecting  layer,  the  layer  tilt  is  included  within  the  iteration 
scheme. 

At  each  stage  of  the  ray-trace  tests  are  made  for  pathological  cases  such  as 
convergence  failures  in  the  allotted  number  of  iterations,  inconsistencies  between 
stages,  extremely  tilted  reflecting  layers,  ducting  of  the  ray  between  layers,  and 
so  forth.  Error  diagnostic  messages  are  printed  prior  to  returning  the  main  pro- 
gram should  TRISL  be  required  to  abort  the  task. 

16.  FUNCTIONAL  DESCRIPTION  C F PROGRAM  OBLFACT 
16.1  Initialization  (Block  1) 

This  part  of  the  program  clears  the  INPUT  file  (NAMELIST  INPUT).  The 
elements  of  the  NAMELIST  are  listed  in  Table  3 including  the  mnemonic  name  used, 
definition,  and  default  value.  A call  to  RHP  clears  the  INPUT  file  of  the  ITS  param- 
eter prediction  coefficients.  A call  to  FACTO  calculates  the  F 2 layer  correction 
factors  from  fQF2  and  h^^Fj  scaled  from  the  vertical  ionogram  and  appearing  in 
the  NAMELIST  input.  Successive  calls  to  PEAK  and  SKIP  with  the  remote  correc- 
tion factors  FR<J  and  MR0,  calculated  from  a previous  run  of  OBLFACT,  check  that 
the  ray  tracing  proceeds  as  expected  (ground-to-ground  mode).  If  not,  the  mission 
is  aborted  and  the  next  task  is  initiated. 


Table  3.  Elements  of  NAMELIST/INPUT/ 


Mnemonic 

Definition 

Unit 

Notes 

SN 

Sunspot  Number 

UNDEF 

YRMD 

Year  & Month,  concatenated 

UNDEF 

DOY 

Day  of  Year 

day 

UNDEF 

ZT 

Universal  Time,  Vercical  Ionogram  (VI) 

hours 

UNDEF 

RI 

Geocentric  Distance,  VI 

km 

6370. 

NL1 

North  Latitude,  VI 

degrees 

WL1 

West  Longitude,  VI 

decrees 
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Table  3.  Elements  of  NA  MR  LIST  /INPUT/  (Cont) 


Mnemonic 

Definition 

Unit 

Notes 

/I  REF 

Universal  T'me,  Oblique  Ionogram  (OI) 

hours 

UNDEF 

RS 

Geocentric  Distance  to  ray  end  point 

km 

5 37  0. 

NLREF 

North  Latitude,  OI  Sounder 

degrees 

WLREF 

West  Longitude,  OI  Sounder 

degrees 

AZ 

Bore  sight  azimuth 

degrees 

E ESTEP 

Elevation  step  for  SKIP  search 

degrees 

1.  0 

FCF2R 

f 1'  from  VI 

muz 

UNDEF 

1 1 MI  NR 

h'  . F_  from  VI 

km 

UNDEF 

min  2 

FP 

Peak  Frequency 

MHz 

DELP 

Peak  Group  Path  Error 

km 

FS 

Skip  Frequency 

MHz 

DELS 

Skip  Group  Path  Error 

km 

DCUK 

Tolerance  in  Group  Path  Deviations 

10.  0 

MR0 

Gradient  Correction  Factor,  M(3000)F2 

4 

F R 0 

Gradient  Correction  Factor, 

MHz 

4 

MODE 

| MODeI  = number  of  hops 

9 

IRSTRT 

Termination  Flag  (STOP  if  0) 

l 

Explanation  of  Terms  Used  in  Table  3 


NOTES  - A numerical  "alue  of  this  column  is  the  default  value  set  in  OBI-FACT . 

UNDEF  - Implies  that  OBLFACT  does  not  establish  a default  value. 

- These  values  are  determined  from  a comparison  between  the  observed 
oblique  ionogram  (001)  and  the  predicted  oblique  ionogram  (POD 
generated  from  a previous  run  of  the  BLOBZ  driving  program.  Refer 
to  "Step  5"  of  he  Operating  Instructions.  Default  values  of  0 are 
assigned  to  F and  DELP. 

4 - For  the  first  run  of  OBLFACT,  these  parameters  are  set  to  0. 

Corrected  values  will  be  part  of  the  output  of  OBLFACT  to  be  used 
as  inputs  of  a new  run  of  OBLFACT  at  the  user  level  of  iteration. 

Refer  to  Steps  4 to  7 of  the  Operating  Instructions. 

? - This  parameter  is  an  element  of  the  argument  list  of  the  ray  tracing 

program  TRISL.  Neither  the  operating  instruction  or  the  program 
flow  chart  illuminate  the  meaning  of  its  sign.  OBLFACT  assigns 
a default  value  of  -1. 
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10.2  KiTuirinriit  of  MR  (Hlork  II) 

This  part  of  the  program  iterates  calls  to  PEAK  in  order  to  refine  the  M(HOOO) 
gradient  correction  factor  Mil.  Since  the  peak  frequency  depends  on  antenna  pattern 
and,  thus  upon  elevation,  it  therefore  depends  primarily  on  the  height  of  the  layer. 
The  Block  II  iteration  maintains  the  fQ F 2 gradient  correction  factor  fixed.  The  value 
of  MR  is  adjusted  initially  by  0.001  and  subsequently  by  a Newton  iterative  technique 
until  the  predicted  peak  group  path  agrees  with  that  from  the  oblique  ionogram  to 
within  the  specified  tolerance  (l)CIIK,  nominally  10  km),  at  which  point  transfer  is 
made  to  Block  III. 

Block  II  is  ultimately  reentered  from  Block  IV  after  refinement  of  the  f F '2 
gradient  correction  factor.  A call  to  PEAK  establishes  whether  further  refinement 
is  required.  If  the  resultant  group  path  is  within  tolerance,  the  task  is  terminated 
and  the  results  printed. 

Should  the  required  information  to  effect  this  refinement  be  not  available,  the 
default  values  of  6 assigned  to  FP  and  DEL.P  automatically  cause  this  part  of  the 
p-ogram  to  be  bypassed.  Thus,  only  the  fQF2  gradient  correction  is  determined; 
the  M(3000)F2  gradient  remains  by  default  the  same  as  that  of  the  parameter  model 
(ITS-7  8),  which,  presumably,  is  the  best  available  in  the  absence  of  evidence  to  the 
contrary. 

16.3  Refirienipnl  <>*'  Fli  (Block  III  and  l\  ) 

This  part  of  the  program  iterates  calls  to  SKIP  in  order  to  refine  the  fQ P' ^ 
gradient  correction  factor  FR.  In  SKIP,  rays  are  traced  at  frequency  FS  from  0° 
to  20"  in  steps  given  by  FLSTFP  (nominally  1°  and  must  not  be  less  than  1/2°).  The 
array  of  group  paths  is  searched  to  locate  the  minimum.  This  is  compared  with 
tin.  minimum  group  path  extracted  from  the  leading  edge  of  the  oblique  ionogram 
at  frequency  FS.  If  the  difference  is  within  tolerance,  the  MR  iteration  is  re- 
entered. If  not,  FR  is  corrected  for  the  nejd  iteration. 

Since  OBI. FACT  is  in  a late  user  stage  rf  simulation  and  a previous  stage 
BLOB  2 produces  the  general  ray-trace  results,  the  user  should  have  limiting  in- 
formation on  the  elevations  between  which  the  minimum  group  path  should  be 
located.  The  built-in  0°  to  20°  may  not  be  sufficient  to  properly  locate  the  mini- 
mum group  path.  The  using  facility  might  consider  modifying  the  program  to  allow 
the  user  to  speci'y  the  elevation  domain  to  be  searched.  This  may  be  accomplished 
as  follows: 

(a)  I. ink  the  initial  elevation,  El, BEG  and  number  of  elevations  to  be  traced, 

N If  I,, VI  to  NA  M El,  1ST  / INPUT/  in  OBI. FACT. 

(b)  Add  FI. BEG  and  NF1.VI  to  COM  MON/ FIDDLE/  in  OBI. FACT.  SKIP, 
and  PEAK. 
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(c)  Modify  the  elevation  DO  loop  in  SKIP  as  follows: 

PRESENT  VERSION  RECOMMENDED  VERSION 

NELM=20,  Dp/ELSTEP  4 1.1)6  (delete) 

EI,=  ELBEG 

DO  10  IEL=1,  NEI.M  DO  10  1EL=1,  NEI.M 

EL=(IEE-1)  !•:  I,  STEP 


10  CONTI N UK 


0 EL=EL  4 E ESTEP 


17.  SI  HIUM  TINK  PEAK  (Fit It,  MUR,  t;i*l\  IGOUAK) 

17.1  Argument  List 

ERR  - fQl"' g correction  quantity  (INPUT) 

M RR  - AKOOOOIF^  correction  quantity  (INPUT) 

GPP  - Predicted  peak  group  path  (OUTPUT) 

IGOHAK  - Error  condition  returned  (OUTPUT) 

I 7.2  Common  Blocks 

/PERT  Cl  linked  to  PERT  F,  the  routine  that  calculates  F0-layer  correction 
factors  from  parameters  scaled  from  the  oblique  ionogram. 

/REMWEN/  linked  from  TRISL.  Contains  reflection  layer  information. 
/FIDDLE/  linked  from  OBI.EACT.  Contains  the  "peak  frequency"  FP  scaled 
from  the  oblique  ionogram,  transmitter  location,  and  take-off  azimuth. 

17..'!  Function 

Computes  the  "predicted  peak  group  path"  corresponding  to  oblique  ionogram 
correction  factors  ERR  and  MRR,  and  "peak  frequency"  FP. 

17.1  Algorithm 

The  input  parameters  ERR  and  MRR  are  stored  in  elements  3 and  5 of  /PERT  C/ 
for  use  by  PERT  E for  calculating  correction  factors  to  fQ F g and  M (3000) F^. 

An  elevation  is  selected  from  the  following  critical  table  according  to  the  peak 
frequency,  element  1 of  /FIDDLE/: 


Table  4.  Elevation  Determination 


FP  less  than 

Elevation 

FP  less  than 

Elevation 

5.  74 

Undefined 

16.  74 

5.  2 

6.  74 

11.  9 

17.  74 

4.  8 

7.  74 

10.  3 

18.  74 

4.  5 

8.  74 

8.  9 

19.  74 

4.  2 

9.  74 

7.7 

20.  74 

4.  0 

10.  74 

6.  7 

21.  74 

3.  8 

11.  74 

5.  9 

23.  74 

3.  7 

12.  74 

5.  2 

25.  74 

3.  6 

13.  74 

4.  6 

28.  74 

3.  5 

14.  74 

4.  1 

30.  74 

3.4 

15.  74 

3.  9 

> 30.  74 

Undefined 

The  calling  program  is  required  to  insure  that  5.  74  < FP  < 30.  74. 

The  ray-tracing  program  TRISL  is  called  to  compute  the  predicted  peak  group 
path  GPP. 

Elements  of  /REMWEN/  are  checked  to  insure  reflection  consistency  occurs  in 
the  F2  layer.  Furthermore,  the  geocentric  distance  of  the  ray  endpoint  is  checked 
that  it  is  less  than  6378.  85  km.  If  these  tests  are  passed,  the  error  condition  flag 
IGOBAK  is  set  to  zero;  otherwise,  an  error  condition  is  noted  by  setting  IGOBAK 
to  unity. 


18.  SUBROUTINE  SKIP  (FRR,  MRR,  GI'S,  10 FI',  IGOBAK) 

18.1  Argument  I.isl 

FRR  - f Fj  correction  quantity  (INPUT) 

MRR  - MOOOOlFg  correction  quantity  (INPUT) 
GPS  - Predicted  minimum  group  path  (OUTPUT) 

IOPT  - Correction  control  parameter  (INPUT) 
IGOBAK  - Error  condition  returned  (OUTPUT) 


18.2  Common  Blocks 

/PERT  C/  linked  to  PERT  F 
/REMWEN/  linked  from  TRISL 

/FIDDLE/  linked  from  calling  program  (OBLFACT) 
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1)1.  3 l''iiiirlion 


Computes  the  predicted  minimum  group  path  corresponding  to  the  oblique 
ionogram  correction  terms  FRR  and  MRR,  and  "skip  frequency"  FS  (element  2 
of  /FIDDLE/). 


1(1.4  Algorithm 

SKIP  transfers  FRR  and  MRR  of  the  argument  list  to  elements  3 and  5 of 
'PERT  Cl  respectively,  to  be  used  in  PERT  F to  calculate  Fg  correction  factors. 

A sequence  of  elevations  between  0°  and  20°  is  selected  for  ray  tracing  via 
calls  to  TRISL.  The  elevation  step  is  contained  in  element  8 of  /FIDDLE/  and  a 
sufficient  number  of  elevations  are  selected  to  include  the  entire  domain  (0°,  20°+). 
Provision  is  made  for  up  to  4 1 ray  traces;  thus,  the  user  must  insure  that  the  eleva- 
tion step  size  is  no  greater  than  0.  5°.  The  elevation  scan  is  terminated  at  the  first 
occurrence  of  a penetration  through  the  ionosphere  if  such  occurs  before  the  normal 
loop  termination. 

The  rays  traced  in  this  elevation  scan  are  partitioned  into  three  sets: 

(a)  Rays  not  reflecting  from  the  F^  layer.  These  are  excluded  from 
further  consideration  in  (b)  and  (c). 

(b)  Rays  with  end  point  at  geocentric  distance  no  greater  than 
6378.  85  km. 

(c)  Rays  with  end  point  at  geocentric  distance  greater  than 
6378.  85  km. 

If  set  (b)  is  empty,  the  error  flag  IGOBAK  is  set  to  unity  and  control  is  re- 
turned to  the  calling  program.  Otherwise,  the  element  of  (b)  with  the  smallest 
group  path  is  determined.  This  least  group  path  is  established  as  the  predicted 
minimum  group  path  GPS  unless  the  following  conditions  obtain: 

(a'  SKIP  had  been  called  with  IOPT  not  zero,  and 

(b)  The  ray  with  next  higher  elevation  belongs  to  set  (c'  and  has  a 
smaller  group  path  length. 

In  the  latter  case,  the  minimum  group  path  returned  as  GPS  is  given  by: 

GPS  = G,  + J(C.  -G,)/(R  -R,)J  (6378.85  - R.  ) , 
b ( c b c b I b 

where  the  G's  and  R's  refer  to  group  paths  and  end  point  geocentric  distance,  and 
the  subscripts  b and  p refer  to  the  above  selected  elements  of  set  (b)  and  (c). 

Should  calls  be  issued  to  SKIP  more  than  100  times  in  any  run  of  OBI.FACT, 
the  error  flag  is  set  to  two  and  the  task  immediately  aborted. 
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19.  HIM  AIIKS  ON  Till.  1 SK  ()!•’  1)01  BI.K  PRECISION  ARITIIMUTIC 


The  version  of  the  program  that  the  present  author  has  studied  has  been 
modified  from  C’lX'  6600  machines  for  use  on  UNIVAC  machines.  The  essence  of 
this  modification  is  to  assign  the  double  precision  attribute  to  every  real  variable. 

In  tlie  opinion  of  the  present  author,  wholesale  use  of  double  nrecision  arithmetic 
is  not  justified,  reduction  notwithstanding. 

In  the  first  place,  double  precision  arithmetic  is  expensive.  Analysis  of  double 
precision  arithmetic  reveals  that  (1)  twice  the  core  storage  is  required  for  the 
variables  involved  in  the  calculation,  (2)  the  coding  requires,  in  general,  three 
times  the  number  of  instructions,  and  (3)  the  execution  time  increases  by  approxi- 
mately a factor  of  four. 

In  the  second  place,  increasing  the  precision  of  a calculation  does  not  increase 
the  accuracy  of  the  computed  result.  For  example,  the  geomagnetic  field  employed 
in  conjunction  with  the  ITS  parameter  prediction  model  is  a sixth  order  spherical 
harmonic  expansion.  The  expansion  coefficients  are  given  to  one  part  per  million 
or  less.  An  eight  significant  digit  machine  in  single  precision  can  calculate  the 
geomagnetic  field  at  a point  to  the  accuracy  of  the  model.  The  ionospheric  param- 
eter model  is  essentially  a one  percent  model;  double  precision  arithmetic  will  not 
improve  this  situation. 

On  the  other  hand,  there  are  situations  where  single  precision  arithmetic  may 
produce  inaccurate  results;  for  example,  in  the  solution  of  plane  triangles  in  cer- 
tain situations. 


Consider  the  situation  depicted  here.  The  distances  R^ 
and  Rg  and  the  angle  0 are  known  and  the  distance  S 
between  Pj  and  P9  is  required.  The  law  of  cosines 
applies  and 

S2  = H1  + H2  " 2H1  R?  cos  9 ■ <30) 

Typically,  the  distances  Rj  and  R^  are  of  the  order 
of  6000  km,  the  arc  Pj  Pj  might  be  60  km  (P  = 0.  01  rad), 
and  the  altitude  difference  Pj  P,;  might  also  be  of  the  order 
of  60  km.  Thus,  take  R9  and  cos  0 to  be  given  exactly  by 

R2  = R (1  a 0.  01), 

cos  0 - 1 - * 10”1  • 
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\\  i tli  single  precision  of  8 significant  digits. 


2 9 -9  -4  -0  -2  -4  -G  -8 

5 = R“  l+n+2'/10  4 10  ± 10  ) - (2+2X10  -10  -10  ± 10  1 

s2  4 r2  ?.y  io'4 *  4 io'6 *  ± io'8  . 

In  this  case,  the  truncation  error  has  resulted  in  reduction  of  precision  to  one  part 
in  104. 

This  loss  of  precision  may  or  may  not  be  tolerable;  in  either  case  it  is  not 
necessary,  for  the  law  of  cos  nes  may  be  reformulated  to  solve  the  triangle  I^P'i 
for  S,  retaining  full  single  precision  results.  The  length  of  the  cord  of  the  angle 
f)  is  given  by 

6 = 2R  j sin  0 . 

Furthermore, 


cos  O'  = - sin-j  0 = -6/21?!  , 

'V  'V  h- 


Then, 

9 9 9 9 

S = h 4 6 4 h6  IH  . (31) 

This  form  of  the  law  of  cosines  will  retain  single  precision  results,  and  should 
be  used  wherever  possible.  To  implement  this  recommendation  in  SUBROUTINE 
TRI.SL,  however,  would  require  considerable  effort  and  analysis;  a totally  new  pro- 
gram may  be  required. 


20.  MODEI.  COMPARISON 


For  purposes  of  comparison,  the  WIMP  program  modelling  subroutines  were 
used  to  generate  ionospheric  models  for  various  combinations  of  parameters  (Sun- 
spot Number,  Season,  Universal  Time)  for  a midlatitude  location.  These  were  com- 

4 

pared  with  similar  models  derived  from  a program  in  use  at  RADC.  These  com- 
parisons arc  shown  in  Figures  3 through  13. 

4.  Rush,  C.  M.  , Miller,  I).,  and  Gibbs,  .1.  ( 1 974)  The  relative  daily  variability  of 

FnF ? and  hmF;j  and  their  implications  for  IIP  radio  propagation.  Radio  Science, 

9:740-756. 

vw 
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Moth  models  use  the  same  formulation  for  f l-'9  the  maximum  plasma  fre- 
quency of  the  I' -layer  (the  ri'S-711  model).  Figures  3 through  5 show  the  Universal 
Time  variation  of  the  height  cf  maximum  electron  density  (h  F,  ' for  three  seasons 
and  two  sunspot  numbers  |\V(XX)  and  H(X.X)  refer  to  the  WIMP  and  KADC  models 
for  SSN  XX].  Similarly,  Figures  6 through  !!  show  the  comparison  of  parabolic 
semithickness  (v^-'^';  WIMP  uses  a semithickness  val  jo  which  is  time  independent. 
'I’lie  WIMP  model  contains  a discontinuity  at  the  lower  .side  of  the  K layer,  which 
results  in  the  most  significant  departures  from  the  RAIX-  model  in  the  vertical  pro- 
file. In  order  to  illustrate  these  departures,  the  difference  between  the  plasma  fre- 
quency was  calculated  by  the  two  models  at  the  r Ititude  of  the  Fj -layer  maximum 
used  in  WIMP.  These  results  are  shown  for  four  seasons  in  Figures  9 through  12, 
where  the  times  of  sunrise  and  sunset  are  also  indicated.  Figure  13  shows  a com- 
parison of  the  two  models  for  given  time  (\  is  the  solar  zenith  angle'.  The  relevant 
ionospheric  model  parameters  are  indicated  in  the  figure. 

The  results  of  ray  tracing  obliquely  through  the  3-0  ionosphere  of  which  f ig- 
ure 13  represents  a single  vertical  profile  are  shown  in  Figure  14.  The  ray-trace 
program  used  was  the  standard  RADC  program,  known  as  ARCON  II,  ITS,  As  ex- 
pected, the  greatest  discrepancies  between  the  oblique  group  path  for  the  WIMP  and 
RADC  models  are  to  be  found  corresponding  :o  tbe  profile  discontinuity  at  the  F ^ 
layer  (compare  Figures  9 through  12).  In  o ■ -d e r to  examine  tbe  differences  in  com- 
puted leading  edge  of  a backscatter  ionogram,  due  to  tbe  different  models,  a one- 
dimensional ionospheric  model  was  generated,  having  the  vertical  profile  given  by 
the  RADC  and  WIMP  models  in  Figure  13,  (that  is,  having  no  horizontal  gradients). 
Leading  edges  were  then  computed  for  simulated  backscatter  ionograms,  using  the 
ARCON  II  ray-trace  program,  and  the  results  shown  in  Figure  15. 
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figure  5.  HA  DC  and  WIMP  Model 
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m 2 


UNIVERSAL  TIME 


figure  (i.  HADC  and  WIMP  Model 
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Figure  9.  RADC  Plasma  Frequency  at 
WIMP  h 1'  Minus  WIMP  f I’.,  Winter 
Solstice 
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Figure  10.  RADC  Plasma  Frequency  at 

WIMP  h F,  Minus  WIMP  f F , Vernal 
„ . m 1 o 1 
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Figure  11.  RADC  Plasma  Frequency  at 

U IMP  h F.  Minus  WIMPf  F , Summer 
c . m 1 o 1 


I'igure  12. 

WIMP  h F 
m 

r.quinox 


RADC  Plasma  Frequency  at 
Minus  WIMP  f F,  Autumnal 
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Figure  13.  RADC'  and  W IMP  Plasma 
Frequency  Profiles 


Figure  14.  A I, CON  I!  IG-.VIHz  Hav  Traces,  HA  DC  ana  WIMP 
Profiles 
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Figure  15.  ARCON  II 
Predicted  Leading  Edges, 
RADC  and  WIMP  Profiles 


21.  COORDINATE  CONVERSION 


The  conversion  from  apparent  radar  range  (that  is,  group  path  (P)  to  true  range 
(R)  can  be  accomplished  if  the  ionospheric  structure  is  known,  by  using  the  proce- 
dures described  earlier.  The  RADC  model  for  1320  UT  on  25  September  1975  was 
used,  together  with  the  ARCON  II  ray-trace  program,  to  compute  (P-Rl  as  a func- 
tion of  P for  frequencies  in  the  range  18  to  26  MHz,  as  shown  in  Figure  16.  The 

locus  of  the  backscatter  leading  edge  is  indicated  by  the  curve  labelled  P . , cor- 

m i n 

responding  to  the  minimum  group  path  at  each  frequency. 

Figure  17  shows  the  comparison  of  the  actual  leading  edge  of  a backscatter 
ionogram  measured  at  1320  UT  (Curve  21  with  the  simulated  leading  edge  (Curve  6). 
Curve  1 shows  a leading  edge  measured  at  1340  UT,  indicating  that  the  ionosphere 
is  moderately  unstable,  since  Curves  1 and  2 are  noticeably  different.  Curve  3 
indicates  a secondary  leading  edge  appeal  ing  in  the  ionograms,  perhaps  a result  of 
azimuthal  gradients  in  the  ionosphere.  The  ionospheric  model  was  then  updated  by 
modifying  it  on  the  basis  of  the  measured  vertical  ionogram  and  allowing  for  a linear 
gradient,  such  that 


N 

e 


BN  ( 1 + A R ) , 
eo  + 


3 9 


where  N is  the  modelled  electron  density  at  any  location. 

Npo  is  the  modelled  electron  density  (from  the  RADC  program) 

(f  F„)  measured 

p _ o 2 

(f  F0)  model  ' 
o i 

A = constant  , 

R+  = boresight  range  (in  radians)  . 

It  should  be  noted  that  the  range  is  here  measured  in  terms  of  the  angle  subtended 
at  the  earth's  center  by  the  arc  on  the  earth's  surface  which  is  the  ground  range. 
Curve  5 in  Figure  17  shows  the  computed  leading  edge  for  A = 2.  0 and  Curve  4 
shows  the  best  fit  (A  = 1.61)  with  the  measured  data.  By  repeating  the  procedure 
illustrated  in  Figure  16,  the  best  estimated  conversion  from  radar  range  to  true 
range  can  be  determined  for  any  operating  frequency.  Although  azimuthal  radar 
information  is  not  available,  it  is  possible  to  ray  trace  at  various  azimuths  to 
establish  a locus  on  the  ground  of  possible  ranges  for  a target  having  any  given  radar 
range. 


Figure  16.  Illustrating  Range  Conversion 
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GROUP  PATH  (km) 


Figure  17.  Backscatter  Leading  Edge 
Comparison 
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Appendix  A 


Program  Listings  for  the  WIMP  Model  Ionosphere 


Program  listings  are  included  for  the  WIMP  model  ionosphere;  each  sub- 
routine is  briefly  described  as  follows; 

1.  R1IP  - Calling  control  program  for  layer  parameters  and 


electron  density. 

2. 

CR  Pi, 2 

- Generates  f F„  and  M(3000)F  . 

O c l 

3. 

MM  DIP 

- Geomagnetic  Field  Model. 

4. 

PERT  F 

- f F,  and  MPOOOlFg  gradient  factors. 

5. 

DFP3  - 

Great  Circle  Range / A zimuth  Calculations 

6. 

DFP  - 

Great  Circle  Range  Calculation. 

7 . 

PTHP  - 

E,  F,,  Y F,,,  and  h Fn  calculations. 

1 m 2 m 2 

8. 

NFROMR  - Electron  Density  Profile. 
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SUBROUTINE  RIIP <RFEC, BNL ATA ,BWL ONA ,PLASC)  RII 

C RII 

C PROGRAMMED  FOR  O.OOOM  BY  COMPUTA TIONS  ,INC . -OCTO BER, 1967  RII 

C RII 

C RII 

C RII 

IMPLICIT  DOUBLE  PRECISION  <A-HvO-Z)  RII 

OOUBLE  PRECISION  H3000*NL AT  A RII 

DOUBLE  PRECISION  MFAC  RII 

LOGICAL  SMOOTH, TILT, NIP,  FAST  RII 

COMMON/RIIPAR/M3000,FCF2,FCF1, FCE , HBE  » HAE , HME,HAF1 * HMF1, HA  F2,  HMF2,  RII 
1SN,ZT,YR,0,HAESET,MFAC,F2FAC,F1FAC,EFAC,H2FAC,SNSQRT,QRAD,CQ,  IDL,I  RII 
2A, SMOOTH  RII 

COMMON/RPERT /RRFEC, IT IP/NIPRIP/  NIP/ T ILTC/TIL T,FA ST  RII 

COMMON/ RIIFHA/IX(4)  RI2 

REAL  IX  RI2 

RI2 

DATA  NL  AT  A,HLONG A RI2 

OATA  SN  ,ZT,YR, 0/13. DO, 13.00, 7509. DO,  2 69. DO/  RI2 

DATA  M3000,  FCF2,FCF1,FCE/ 3.0500, 8.D0, 4.  OBD  0,2.7900/  RI2 

DATA  HBE  ,HHE  ,HMFl,HMF2/90  . DO  , 110  . DO , 18  0. DO,  20  0.  DO/  RI2 

OATA  HAE,HAFl,HAF2/20.D0,70.O0,108.D0/  RI2 

OATA  MF AC , F2FAC , F1F AC ,EF AC , H2FAC/5*1. DO/  RI2 

DATA  HAESET/15.D0/  RI2 

RI2 

NLATA=BNL ATA* 57. 2957795130023200  RII 

HLONGA=BHLONA  *57.2957795130823200  RII 

RRFEC=RFEC  RII 

NL ATA=OMOO( NLAT A, 90 .DO)  RII 

HLONGA=DMOD(WLONGA,360.O0)4180.D0-DSIGN(180.D0,HLONGA)  RII 

1004  CALL  CRPL2( NL  AT  A, MLONGA)  RII 

M3000=M3000*MFAC  RI2 

FCF2  = FCF2*  F2FAC  RI2 

C PRFOICTEO  .BOTTOM  OF  E LAYER  RII 

1005  BA= 10 0 . DO  RII 

AA=BA4HAESET  RII 

IF ( IDL* IA.EO . 0 ) GO  TO  1006  RII 

CALL  DFP(NLATA,WLONGA,76.DO,102.DO,DIST>  RII 

BA=100.D0-4.D0*DHAX1( (1. DO-DABS  CD 1ST- 3.D3)/3.D3)**2,0.0D0)  *DMAX1C  RII 

1 OCOS((DABS(15.D0*DABS(ZT-12.O0)-MLONGA>-90.D0)/(2.O0*  RII 

2 57.29577951308232D0) ) , u . 0D0)  RII 

C PREDICTED  TRUE  HEIGHT  PROFILE  RII 

1006  CALL  PTHP (NLATA,HLONGA,AA,BA)  RII 

C GET  ELECTRON  DENSITY  AT  REF.  POINT  RII 

C GET  N FROM  R RII 

IF ( RFEC .LE. 6370. DO)  RETURN  RII 

CALL  NFROMR(RFEC,PLASD)  RII 

RETURN  RII 

ENO  RII 
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SUBROUTINE  CRPL2( NLAT, HLONG)  CRP 

C ...TO  FINO  CRITICAL  FREQUENCY  AS  A FUNCTION  OF  THETA  AND  PHI  CRP 

C FROM  NEW  N8S  IONISPHERIC  PREDICTIONS  FOR  A PARTICULAR  TIME...  CRP 

IMPLICIT  OOUBLE  PRECISION  (A-H,0-Z>  CRP 

DOUBLE  PRECISION  NLmT ,M3 000 , MFAC  CRP 

DIMENSION  SINTO(U)  ,COSTO(il>  ,G(76) ,GAMMA(2) ,DKT0<76,2)  CRP 

DIMENSION  CT(8),ST(8>  CRP 

COMMON/RIIPAR/M3000,FCF20,DUN(9> ,SN, ZT , YR , DAY ,HAESET ,MFAC , FFAC,  CRP 
1DUM1  ( 6)  ,10(3)  CRP 

*/LABCRP/NHARM(2),KI(2),KII(2>  , Kill ( 2) , KIV(2> , KV (2) , KVI (2  ) , KV II ( 2 ) , CRP 
*KVIII(2) , KIX(2) ,0(13,76,2)  CRP 

DATA  RAD/. 01745329251 99433D0/  CRP 

G(l)=1.  00  0 CRP 

CALL  MMOIP(NLAT,WLONG,SSMX)  CRP 

THERAO= (360.O0-WLONG)*RAO  CRP 

COSLAM=DCOS(NLAT*RAO)  CRP 

CT(1) =COSLAM*OCOS (THERAO)  CRP 

ST( 1) =C  OSL AM  *DSIN (THERAD)  CRP 

00  63047  IJK=2 , 8 CRP 

CT  ( IJK)  =CT  ( 1 ) *CT(  IJK-l)-ST(l)  *ST(IJK-1)  CRP 

63047  ST(IJK)  = CT(1)*ST  ( IJK-1) *ST( 1) *CT (IJK-1)  CRP 

IF(DABS(ZT-ZTP) .LT. 1.0-3.  AND.  DABS ( YR-Y RP) . LT . 1 . D-l) GO  TO  20  CRP 

ZTP=ZT  CRP 

YRP=YR  CRP 

T 0= 15. 0 DO*  ZT-180 • 0 00  CRP 

T0RAD=T  0*RA0  CRP 

ARG=0 .0  DO  CRP 

NHARML=MA  XO (NHARM ( 1) , NH ARM ( 2)  ) CRP 

DO  1 J= 1 , NHA  RML  CRP 

ARG=T  ORAD  + ARG  CRP 

SINTO(J)=DSIN(ARG)  CRP 

1 COS TO (J)=OCOS(ARG)  CRP 

DO  2 L= 1, 2 CRP 

NHARML=  NHARM (L)  CRP 

K9=KIX(L)+1  CRP 

DO  2 K=  1,K9  CRP 

DKTO(K,L)=D(l,K,L)  CRP 

00  3 J=1,NHARML  CRP 

3 OKTO( K, L) =OKTO( K, L) *0 (2* J+l , K,L)  *COSTO ( J) *0 ( 2*J,K,L) *S INTO ( J ) CRP 

2 CONTINUE  CRP 

20  DO  100  L=l, 2 CRP 

NHARNL= NHARM (L)  CRP 

K1=KI  (L)  *1  CRP 

K2=KII(L)*1  CRP 

K3=KIII(L)*1  CRP 

K4=KIV ( L) * 1 CRP 

K5=KV(L)*i  CRP 

K6=KV1(L)  *1  CRP 

K7=KVII  (L  ) *1  CRP 

K8=KVIII(L) *1  CRP 

K9=KIX(L) *1  CRP 

KIP 1=K1 *1  CRP 

K1P2=K1 *2  CRP 

K1P3=K1*3  CRP 

K2P 1=K2*1  CRP 

K2P2=K2  *2  CRP 

K2P3=K2+3  CRP 
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K 3P  1=K  3 + 1 CRP 

K3P2=K3+2  CRP 

K3P3=K3«-3  CRP 

K4P1  = K4  +1  CRP 

K4P2=K4+2  CRP 

K4P 3=  K4  *3  CRP 

K5Pi=K5  +1  CRP 

K6P1=K6*1  CRP 

K7P 1=  K7* 1 CRP 

K8P1=K8  +1  CRP 

00  4 K=2,K1  CRP 

4 G<K)=G(K-i)*SSMX  CRP 

IE(  Kl.EO.  K2)  GO  TO  9 CRP 

G (K IP  1) =CT  (1 ) CRP 

G(K1P2)  = ST(  1)  CRP 

00  5 K=K1P3,K2,2  CRP 

G ( K)  =G<  K-2)  *SS*1  X CRP 

5 G(K*1)=G(K-1)*SSMX  CRP 

IF  ( K2.E  Q.  K3)  GO  TO  9 CRP 

G(K2P1)  =CT  (2)  CRP 

G ( K2P2) =ST ( 21  CRP 

00  6 K=K2P3,K3,2  CRP 

G(K)=G(K-2)*SSMX  CRP 

6 G (K*l)  = G( K-i) *SSMX  CRP 

IF(K3.EQ.K41  GO  TO  9 CRP 

G(K3P1)  — C T ( 3 ) CRP 

G (K  3P2)  =ST  (3 ) CRP 

00  7 K=K3P3,  K4,  2 CRP 

G(K)=G(K-2)*SSMX  CRP 

7 G(KU)  = G(K-11  *SSMX  CRP 

I F(  K4« E 0.  K5>  GO  TO  9 CRP 

G (K 4P  1)  =CT  ( 4 ) CRP 

G ( K4P2)  =ST ( 4)  CRP 

00  8 K=K4P3,K5,2  CRP 

G(K)=G<  K-2)  *SSMX  CRP 

8 G ( K + l)  =G  { K-l ) *SSMX  CRP 

IF(K5.EQ.K6)  GO  TO  9 CRP 

G ( K5P1)  =CT  <5  ) CRP 

G (K6)  = ST  (5  ) CRP 

IF(K6.EQ.  K7)  GO  TO  9 CRP 

G (K6P1)  =C  T (6  ) CRP 

G ( K7) =ST( 6 ) CRP 

IF (K7.E  Q. K8)  GO  TO  9 CRP 

G(K7P1) =CT  (7)  CRP 

G(K8)=ST(7)  CRP 

IF(K8.EQ.K9)  GO  TO  9 CRP 

G ( K8P 1) =CT  ( 8 ) CRP 

G (X9)  =5  T ( 8 ) CRP 

9 GAMMA(L)=0.000  CRP 

00  10  K=1 ,K9  CRP 

10  GANMA(L)=GAMMA(L) +OKTO(K,L) *G(K)  CRP 

100  CONTINUE  CRP 

FCF20=GAMMA( 1) *0ABS(FFAC)  CRP 

M3000  = GAMHA(2)*OABS  (MFAC)  CRP 

IFtFFAC.LT.O .OOOJCALL  PERTF ( 7 Y, NL A T,  HL ONG , FCF 20>  CRP 

IF(MFAC.LT.O.OOO)CALL  PERTH ( ZT , NLA T , WLONG, H3 0 00 ) CRP 

RETURN  CRP 
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SUBROUTINE  PTHP ( NL AT ,HLONG, A, B)  PTH 

C PREDICTED  TRUE  HEIGHT  PROFILE  PTH 

C PROGRAMMED  FOR  O.ODOM  BY  I.STUHLER  PTH 

IMPLICIT  OOUBLE  PRECISION  <A-H,0-Z>  PTH 

OOUBLE  PRECISION  NL AT ,M3 000 , H, J, KK  PTH 

LOGICAL  SMOOTH  PTH 

COMHON/FACCOM/FM,ZZ,KK,P3000  PTH 

COMMON/  R IIP  AR/ M3  0 00  * FCF2  Q,FCF10,FCEO,HBEP»HPE»HHEP*',rFl»  HHF  IP  ,HPF  2 PTH 
1.HMF2P,  SN,ZT,YR,0,OUM(3)  ,F< FAC,EFAC, H2FAC,  PTH 

2 SNSQRT,QRAD,CQ, 10(2), SMOOTH  PTH 

OATA  RAO/57. 2957795130823200/  PTH 

C SOLVE  FOR  ZENITH  CORR.  ANGLE  FOR  DAY  OF  YEAR  PTH 

M = 23.  45O0*OSIN(  180.00 *<0-82.500) / ( 1 8 2 . D0*RA 0) ) PTH 

C GET  F 30 00  PROPAGATION  FREQ  FOR  F2  LAYER  PTH 

F3000=FCF20*M3000  PTH 

C CHECK  FOR  INCORRECT  OAY  (NEG)  PTH 

IFIO.LT. 0.000)  GO  TO  2000  PTH 

C IS  DAY  PAST  SPRING  PTH 

IFCO.GE. 162.500)  GO  TO  1000  PTH 

X=0+365 ,500  PTH 

Y = l.  ODD  PTH 

GO  TO  1001  PTH 

C IS  OAY  PAST  FALL  PTH 

1000  IF(0.GE.344.  500  ) GO  TO  1002  PTH 

X=0  PTH 

Y=-i. GOO  PTH 

GO  TO  1001  PTH 

C CHECk  FOR  INCORRECT  OAY  (TOO  LARGE)  PTH 

1002  IF(D.GE.366. 500)  GO  TO  2001  PTH 

X=D  PTH 

Y=1 . 00  PTH 

C SOLVE  FOR  SOLAR  ZENITH  ANGLE  IN  RAO  PTH 

10  01  J = Y*. 023500+ (64.0  0/  (4 .0 0+0 ABS ( ( X- 344,  5D0)  /45.00)  **3) -240 .DO/ (15.00  PTH 

1 +OABS(X-344.500) ))  PTH 

IF(ZT.LT. 12.00) ZT=-ZT  PTH 

Z=DABS(  15.00 *OA  BS ( ZT-12 • OO ) -J+15 .OQ-WLONG)  PTH 

IFCZ.GE.  180.  00)  Z=  360.  DO-Z  PTH 

ZSQ=Z*Z  PTH 

QSQ=ZSQ*(0C0S (0.  600 + ( NL AT+M) / RAD)  + *2)  + ( NLAT-M) +*2  P TH 

1005  Q=  DSQRT  ( QSQ)  PTH 

QRAQ=Q/ RAO  PTH 

CQ=OCOS(.  932O0+QRAO)  PTH 

C FUNCTION  FOR  E-LAYER  NIGHT  TIME  FORMATION  PTH 

CCQ=( (CQ+1.D0) *.500)**2  PTH 

SNSQRT=OSQRT (1.D0+.004D0+SN)  PTH 

IF(CQ.LE.O.OO)  GO  TO  1101  PTH 

C SOLVE  FOR  PREO  E-LAYER  CRITICAL  FREQ  PTH 

C 10.0489=3.  17*3.  17,  . 56 25= . 75*.  75  PTH 

FCEO=SNSQRT  + DSQRT  (10,  048 9D0*CQ**( 2.  00/3.77500)  + . 562 5D0*CCQ*+ 2)  PTH 

IF(FCEO.GE.FCF20) FCEO=. 75O0*SNSQRT*CCQ  PTH 

GO  TO  1102  PTH 

1101  FCEO=.75DO*SNSQRT*CCQ  PTH 

C SOLVE  FOR  PREO  Fl-LAYER  CIRTICAL  FREQ  PTH 

1102  FCF10  = 1.2600  *FCEO+0  .500  PTH 

IF(FCF10.LT.  FCF20*,  99D0  ) GO  TO  9982  PTH 

FCF 10=. 9900*FCF2Q  PTH 

FCEO=OMAXl ((FC FI 0-. 500)/ 1.2600,1.0- 5)  PTH 
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C SOLVE  FOR  REF.  FREQ.  NO.  1 PTH 

9982  FK1=FCF20-(FCF20-FCF10)/4,DO  PTH 

C CALCULATE  PREO.  SEHITHICKNESS  OF  F2  LAYER  PTH 

HPF2=(0.48900*SN+56.000)*H2FAC  PTH 

C SOVE  FOR  REF  FREQ.  NO.  2 PTH 

FK2  = FCF20  Ml  • O0-(  OSQRT(HPF2)  /10 0.000)  ) PTH 

C SET  FM  TO  GREATER  OF  TWO  REF.  FREQ.  PTH 

IF(  FK1.GT • FK2)  GO  TO  1003  PTH 

FN=FK  2 PTH 

GO  TO  1004  PTH 

1003  FM=FK1  PTH 

C SOLVE  FOR  VIRTUAL  HT.  OF  F3000  PTH 

1004  XX=FH/F30Q0  PTH 

H3000=3100.0Q*XX»+2+70.000  PTH 

FC=4. OO  +FM-3 .D0*FCF10  PTH 

P3000=( (FC+3. 78DO*FCEO+1.5DO)/(1.26DO*FCEO+.5DO) ) *OLOG < < FC+8 . 8200  PTH 

1 *FCEO+3.5DO)/(FC-1.26DO*FCEO-0.5DO))  PTH 

ZZ=(A-B)*  CFC+3.78DO*FCEO+ 1.500)/ (8. DO*FCEO ) +OLOG ( ( FC+7 .78D0  PTH 

* * FCEO+1  • 500  ) / ( FC-.  2 200  +FCEO  + 1.5DO)  ) PTH 

FF=FM/FCF20  PTH 

KK=HPF2*FF*0L0S( (l.OO+FF) /(1.D0-FF) )/ 2.000  PTH 

C CALCULATE  PREO.  HT.  OF  F2  LAYER  MAX.  PTH 

HMF2P=8 .DO*(H3000-ZZ-B-KK)/P3000+HPF2+A  PTH 

C CALC.  HTS.  OF  OTHER  LAYERS  PTH 

HHF1P=HMF2P-HPF2  PTH 

HMEP=  A PTH 

HBEP=B  PTH 

HPF1=HMF1P-HNEP  PTH 

HPE= A-8  PTH 

ZT  = OABS ( ZT)  PTH 

FCEO=OMINi(FCEO*EFAC, <.99DO*FCF20-.5DO)/1.26DO)  PTH 

IF(FIFAC)  1014, 1014, 1015  PTH 

1014  FCF10=1.2600*FCEO+.500  PTH 

GO  TO  1016  PTH 

1015  FCF10=FCF1Q*F1FAC  PTH 

1016  IF(SMOOTH)FCF10=FCF10/. 86602540378443900  PTH 

FCF10=OMIN1( FCF10,  ,9900*FCF  20)  PTH 

RETURN  PTH 

2000  WRITE (6,4)0  PTH 

4 FORMAT ( 4H  0 =013.5, 12H  IS  NEGATIVE)  PTH 

STOP  PTH 

2001  WRITE (6,5)0  PTH 

5 FORMAT ( 4H  D =D13.5,22H  IS  GREATER  THAN  366.5)  PTH 

STOP  PTH 

ENO  PTH 
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SUBROUTINE  NFROMRtRFEC ,PLAS0)  NFR 

C NFR 

IMPLICIT  DOUBLE  PRECISION  <A-H,0-Z>  NFR 

OOUBLE  PRECISION  K,J,M3000  NFR 

LOGICAL  SMOOTH  NFR 

COHMON/RIIPAR/M3QOO  , FCF2 , FCF1, FCE , HBE , KAE , HME , HA  FI , HMF  1,HAF2,HMF2,  NFR 
lDUM(lO)  .SNSQRT, QR  AO  , CO,  ID  , I A , SMOOTH  NFR 

C HEIGHT  ABOVE  GROUNO  OF  REF.  PT.  NFR 

OH=RFEC-6370.D0  NFR 

C NO  - IS  PCINT  OF  INTEREST  ABOVE  E LAYER  MAX  NFR 

IF (OH.GT.HME)  GO  TO  201  NFR 

IF(IO.EQ.l)  GO  TO  2001  NFR 

C IS  PT.  OF  INTEREST  BELOW  E LAYER  NFR 

IF((HME-OH)  .GE.HAE)  GO  TO  300  NFR 

C NO  - SOLVE  FOR  E LAYER  PLASMA  DENSITY  NFR 

F02=FCE*FC£*  <1. D0-< <HME-OH> /HAE)  **2>  NFR 

GO  TO  202  NFR 

C E LAYER  WITH  O LAYER  TAIL  NFR 

2001  ZP=2.D0*(HME-OH)/HAE  NFR 

K= 1. 5 00*0 SORT  Cl. 700/(1. 7D0  + ZP))  NFR 

J=1.D0+ (ZP-1.D0) *OABS(ZP-l.D0)/10.D0  NFR 

J=J/<1.00-. 100/(1. D0+ZP)**2>  NFR 

ZPP=2.00  NFR 

IF(ZP.GT.O.DO) ZPP=2.00-OLOG(1.DO+1.72DO*ZP+*K)  NFR 

IF  (ZPP. LE . 0 . DO ) GO  TO  300  NFR 

ZPP=2.D0-ZPP**J  NFR 

F02=FCE*FCE* ( 1 . DO -ZPP+ZPP/4 . DO)  NFR 

GO  TO  202  NFR 

C SOLVE  FOR  REF.  HTS • 1 THRU  3 NFR 

201  HS2=HMF2-HAF2*OSQRT  ( 1 .00- ( FCF1/FCF2)  **2)  NFR 

HS1  = HMF1-HAF1  + 0SQRT  (1  .0 0- (FCE/FCF 1)** 2)  NFR 

HT2=(HS2+HMFl)/2.O0  NFR 

C IS  PT.  OF  INTEREST  ABOVE  REF  PT  3 NFR 

IF(OH.GT.HTE)  GO  TO  203  NFR 

C NO-  SOLVE  FOR  FI  LAYER  PLASMA  DENSITY  NFR 

BIGK1=(HS1-HME+HS2-HMF1)/(2.D0*(FCF1-FCE) ) NFR 

BIGK3=HAF1/FCF1  NFR 

BIGK2=HMF1- (FCF1*  (HS1-HME)  + FCE* (HS2-HMF1)  ) / ( 2.D0*(FCFl-rCE) ) NFR 

BIGA=BIGK1**2+BIGK3**2  NFR 

BIGB=2.00+8IGK1*  (BIGK2-OH)  NFR 

BIGC=OH**2-2.DO*OH*BIGK2+BIGK2**2-HA  F1+  + 2 NFR 

F 0= (-BIGB+OSQRT ( BIGB* *2-4. D0*BIG A* BIGC) ) / ( 2 . O0*B IG A ) NFR 

IF(FO.LT.FCE) F0=FCE  NFR 

F02=F  0*F0  NFR 

GO  TO  202  NFR 

C SOLVE  FOR  OEMS  IT  Y ABOVE  F2  LAYER  MAX  NFR 

211  OHH= (OH-HMF2) /OH  NFR 

Z=(OH-HMF2)/HAF2*2.DO  NFR 

F02=FCF2**2*0EXP  ( 1 .00 -Z-DEXP (-Z ) ) + ( 0 . 5D0  + FCF2/10 . 0 DO ) *OHH*HMF2/OH  NFR 
F02  =DMIN1  ( F02 i FC F2* *2 ) NFR 

GO  TO  202  NFR 

C IS  POINT  OF  INTEREST  ABOVE  F2-LAYER  MAX  NFR 

203  IFIOH.GT. HMF 2)  GO  TO  211  NFR 

C NO-SOLVE  FOR  F2-LAYER  PLASMA  DENSITY  NFR 

BIGK1=(HS2-HMF1) / (2 . DO* ( FCF2-FCF1) ) NFR 

BIGK3=HAF2/FCF2  NFR 

8IGK2=HMF2-FCF2*8IGK1  NFR 
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8IGA  = BIGK1**2*BIGK3**2  NFR 

BIGB=2.00*BIGK1MBIGK2-OH)  NFR 

8IGC=0H**2-2  .DO*OH*BIGK2*BIGK2**2-HAF2**2  NFR 

FO=  (-8IG8+OSQRT(8IG8**2-4.DO*BIGA*8IGC))/(2.DO*BIGA>  NFR 

F02=F0*F0  NFR 

202  PL A SO=F  02* 12400. 00  NFR 

IF(IO.EQ.l)  GO  TO  301  NFR 

RETURN  NFR 

300  PL  ASD=0 .DO  NFR 

IF(IO.EQ.O)RETURN  NFR 

301  CQ=OCOS(Q RAD*. 8500)  NFR 

IFICQ.GT. 0.000)  GO  TO  302  NFR 

CQ=0. DO  NFR 

RETURN  NFR 

302  CO=CQ** (0.2500)  NFR 

30  3 PL ASD=PLASD*SNSQRT*  SNSQRT*CQ*1 . 02/(1. DO*81.DO*(  (65 • OO-OH) /HAE) **  4)  NFR 

RETURN  NFR 

ENO  NFR 
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SUBROUTINE  OFP ( NLA TA  .WLONGA  ,NLA TB , HLONGB,  OIST)  OFP 

IMPLICIT  OOUBLE  PRECISION  CA-H.O-Z)  DFP 

OOUBLE  PRECISION  NLATA.NLATB  DFP 

DATA  RAO/57. 29577951308232D0/  DFP 

DACOS (DUMMY  )=DAT AN  2 (OSQRT< 1.00 0-DUMHY+DUMMY) .DUMMY)  DFP 

C SHIFT  INPUT  LAT.  BY  ONE  QUA  ORANT  OFP 

PPA=90.  ODO-NLAT  A OFP 

PPB=90.000-NLATB  DFP 

GGH=HLONGA  DFP 

GGK=HLONGB  OFP 

C FINO  SMALLEST  ANGULAR  DIFF . BETWEEN  A AND  B OFP 

HHGK=GGH- GGK  DFP 

ABHHGK=  DABS<  HHG  K)  DFP 

HHK=ABHHGK  DFP 

IFCHHK.GT. 180,00)  HH K= 360 .D 0-HHK  DFP 

C SOLVE  FOR  BEARING  FROM  A TO  B OFP 

RAAPB=HHK/ RAO  DFP 

C SA AP B=OCOS ( RAAP  B)  OFP 

RPPA= PP A/RAO  DFP 

S INPPA=  OSIN( RPPA)  DFP 

COSPPA=  OCOS (RPPA ) DFP 

RPPB=PPB/RAO  DFP 

SINPPB= OSIN (RPPB)  OFP 

COSPPB=OCOS( RPPB)  OFP 

COS  AAB=COSPPA*COSPPB+SINPPA*SINPPB*CSAAPB  DFP 

RAAB=OACOS(COSAAB)  DFP 

C SOLVE  FOR  GREAT  CIRCLE  OISTANCE  BETWEEN  A AND  B DFP 

DIST=637D.0D0*RAAB  DFP 

RETURN  DFP 

END  DFP 
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ROUTINE  HMOIP 


74/74  OPT  = l 


FT  N 4.5*414 


04/27/ 


SUBROUTINE  MMDIP <NLAT,WLONG,SSMX>  DMM 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z)  DMH 

COMPUTE  MAGNETIC  FIELD  COMPONENTS  MAGNETIC  DIP  AND  MODIFIED  MAGNETIC  DMM 
DIP  DMM 

OOUBLE  PRECISION  NL AT  DMM 

DIMENSION  P<7,7)  ,DP<7,7)  , C P (7  ) , A OR  (7  > , SP  (7  ) DMM 

DIMENSION  CT(7,7),H(7,7),G<7,7)  DMM 

DATA  P,0P,SP/1.  000,104*0.  BOO/, CP/i.  000,6*0.000/  DMM 

DATA  RO/57.  29577  95  13  082  3 2D  0 / ,,  HC/6371 . 20  3/  OMM 

DATA  CT/2*0. DO,  .3333333300  , .2666666700, .2571428600  , DMM 

1 .2539682500 2525252500,3*0. CO,. 200,, 2285714200,. 238 095 2 3 DO,  OMM 

2 . 2424242400,4*0.00, . 1428571400, . 190  4761900  ,.  2121212100,  DMM 

3 5*0.00  ,.1111111100  , .161 61 616 DO, 6*  0.00, .0909 0909D0,  14*  0.00/  DMM 

DATA  G/0. 00, .30411200  , . 02403500  , -.03151800, -.  041794D0,  DMM 

1 .016256D0, -.01952300, 0.00  ,.  02147400  , -.051253D0,  OMM 

2 .06213000, -.04529800,-. 03440700, -.00485300, 2*0.00,  OMM 

3 -.013381D0,-. 02489800, -.02179500, -.01944700, .00321200,  DMM 

4 3*  0.  00  ,-.00  649600  , .00700800, -.0G0608D0,.  02141300,  OMM 

5 4*0. DO, -.  00204400,  .00277500, . 001051D0, 5*0. DO, . 00069700,  DMM 

6 .00022700,6*0.00,. 00111500/  OMM 

DATA  H/8*0. DO, -.  05798900, . 03312400,.  01487000, -.  0118  25D0,  DMM 

1 -.  00079600, -.  00575800, 2*0  . DO ,- . 00 1 5 79D0 , - .0  04 075  DO , OMM 

2 . 01000600, -.00200, -.  00873500,  3*0. DO,. 00021DO,. 0004300  , OMM 

3 . 00459700,-. 00340600, 4*0. DO, . 001385DO,. 00242100,  DMM 

4 -.00011800,5*0.00,-.  00121800, -.00111600, 6*0. DO,  DMM 

5 -.00032500/  OMM 

P 1 = NL  AT  DMM 

P2-360.0D0-( WLQNG)  DMM 

IF(P1-89.9D0)2,4 ,1  OMM 

1 P 1=  89. 9D0  DMM 

P 2=  0. 00  DMM 

GO  TO  4 DMM 

2 IFIP1+89. 900) 3,4,4  DMM 

3 Pl=-89. 9D0  DMM 

P2=O.DO  DMM 

4 PHI=P2/RD  OMM 

AR=HC/(HC*3.  D5)  DMM 

C=0SIN(P1/R0)  DMM 

S=OSQRT (l.-C**2)  DMM 

SP ( 2) =DSIN (PHI)  DMM 

CP  ( 2)  =OCOS  (P  HI ) OMM 

AOR(l)=AR**2  OMM 

AOR(2)=AR**3  DMM 

DO  5 M=  3, 7 DMM 

SP(M)  = SP(2)*CP(M-1)+CP(2)*SP(H-1)  DMM 

CP(M)  = CP(2)*CP(H-1)-SP(2)  *SP(M-1)  DMM 

5 AOR(M) =AR*AOR(M-l)  DMM 

BV=0. OO  DMM 

BN=0 • OO  DMM 

BPHI=0. OO  OMM 

OO  6 N=  2, 7 DMM 

FN=  OBLE  (FLOAT(N) ) DMM 

SUMP=  0 • DO  DMM 

SUMT  =0. OO  DMM 

SUHP=0.00  DMM 

OO  7 M=i,  N DMM 

IF(N-M)  8,9,8  DMM 
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74/74  OPT  = 1 


ETN  4.5+414 


04/27/ 


9 P(N,N)=S*P(N-1,N-1>  DHM 

OP(N,N)=S*OP<N-l ,N-1)  +C*  P(N-l,N-i)  OMH 

GO  TO  10  DHH 

8 P(N,M)=C*P(M-i,M)  -CT(N,M>  *P<N-2,M>  OHM 

OP(N,H)=C*OPCN-l,H»-S*PCN-l,H)-CT(N, M)*DP(N-2,M>  DMM 

10  FM=08LE(FL0AT(M-1) ) DMM 

TS=G(N,H)*CP(M) +H(N,M)*5P(H)  DMM 

SUMR=SUMP+P(N,M) *TS  OMH 

SUM  T=  SUHT  +DP (N ,M ) +TS  DMM 

7 SUHP=SUMP  + FM*P(  N,  Ml  *< -G(N,MI  *SP(M>  +H<N,H>  *CP  (M)  ) OMM 

9V=  8V  + A0R (N1  *FN*SUMR  DMM 

BN=BN-AOR(NI  *SUMT  DMM 

6 8PHI=BPHI-A0R(N)*SUMP  DMM 

COSLAM=DCOS ((NLATI/RO)  DMM 

C0M1=-BV  DMM 

COM  2= BN  DMM 

COM  3=-BPHI/S  DMM 

TMP=COM2** 2+C0M3**2  OMM 

BIGI  = 0ATAN(-C0M1/DS0RT(TMP>  I DMM 

SSHX=(BIG I/O  SORT (BIGI**2+COSLAMI ) DMM 

NEXT  CARO  OELETEO  OSINCE  SMX  NOT  USED  BY  CRPL2  OMM 

SMX=OASIN(SSMX)  DMM 

RETURN  DMM 

END  DMM 
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IROUTINE  PE.’-F  74/74  OPT=l  FTN  4.54-414  04/27/ 

SUBROUTINE  PERTFCT,  NL,HL,  FCF20)  PER 

IMPLICIT  DOUBLE  PRECISION  CA-H,0-Z>  PER 

DOUBLE  PRECISION  NLREF, HR , HZ, MO, NL  PEP 

COMHON/PERTC/NLREF,  WL  REF,  FR,FO,  MR,  MD,FZ,MZ,  OUT,  DW,X,  GAMMA,  AZM,  PER 

1 OF  AC , DNAC , A ZREF  PER 

RTD=57. 2957795130823200  PER 

CALL  DFP3(NLREF,WLREF,NL,WL,DIST,AZ)  PER 

AZ= AZ-AZREF  PER 

STHINV=1 .0D0/OABS{DSINICAZM-(90.D0+GAHMA) ) /RTO)  ) PER 

y=OUT*36liO.OO*X*STHINV  PER 

H=OH*STHINV  PER 

SO=OSIN  ( (OIST+Y)  /H*  36  0,  00 /RTO)  PER 

FCF20=FCF20M1.00+FZ*SD)  * ( l.DO*FR*DIST/FO)  M 1 .00*AZ*DFAC>  PER 

FG'/20=AMAX1<  ,51D0,FCF20>  PER 

RETURN  PER 

ENTRY  PERT H(T  , NL ,W  L , FCF20  ) PER 

FCF20=FCF2Q*  (i.OO+MZ*SD> ♦ { 1.  DO  + KR*DIST/MO)  *(  1 .D0+AZ*DMAC)  PER 

FCF20=AMAXl(i,D0,FCF2  0>  PER 

RETURN  PER 

END  PER 
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ROUTINE  FACTO 


74/74  OPT  = 1 


FTN  4.5+414 


04/27/ 


SUBROUTINE  F ACTO ( NL  REF, WLREF, FCF2R, H MINR, ZTREF)  FAC 

IMPLICIT  00U9LF.  PRECI  CION  ( A-H,0-Z>  FAC 

DOUBLE  PRECISION  M300  OP,  M30  0 OR , MF  AC,  MM  , NL  REF , K1  FAC 

COMMON/ FACCOM/F’i,ZZ,KK,P30  0 0 FAD 

COMMON/RIIPAR/:13  000P,FCF2,FCF1,FCE,  HBE,  HAE,  HME,  HAF1 , HMF1  ,HPF2,  FAC 

1 HMF2,SNS<5)  , MFAC,  F2FAC,  0UM<6>  , IDUM<  3)  FAC 

DATA  RTO/57.  2957795130823200  / FAC 

F2FAC=l.O0  FAC 

MFAC=i. 00  FAC 

IF(FCF2R.EQ.O.OO)  RETURN  FAC 

ZTS=SNS<2>  FAC 

SNS  (2)  = ZTREF  FAC 

CALL  CRPL2(NLREF,HLREF)  FAC 

F2  F AC=D A8S  CFC  F2R/FC  F2 ) FAC 

IF  (HMINR.EQ.O.OO)  GO  TO  10  FAC 

CALL  RIIP(  6370.D0,NLREF/RTD,WLREF/RTD,PLASD>  FAC 

Pl=  <FCF2+3.78D0*FCE+ 1.500)  /FCF i*DLO G ( < FCF2  + 8 . 820  0*FCE+ 3.  500)  / FAC 

1 ( FCF2-FCF1)  ) FAC 

MM=HAE/fl.OO  FAC 

Zl=MM*(FCF2+3.78O0*FCE+i.5DO)/FCE*DLOG( (FCF2  + 7. 780  0 *FCE+ 1 . 500 ) / FAC 
1 ( FCF2- .2200 +FCE  + 1 • 500) ) FAC 

FFM=FCFl+fFCF2-FCFl)/4.D0  FAC 

FF 1=FFM/FCF2  FAC 

K1=50.DO*FF1*OLOG((1.00+FF1>/ (1.D0-FF1) ) FAC 

H30  00=P30  0 0/P1*(OABS(HKINR)  -Zl-HBE-HPF 2*K1/1 0 0 .D  0 ) + Z Z+ HBE+ KK  FAC 

M 30  00  R=  (FM/FCF2  )*  OSORT (3100.00/CH3000-70.DO)  ) FAC 

MFAC=DSIGN(M3000R/M3000P,HMINR)  FAC 

10  F2FAC=0SIGN(F2FAC,FCF2R)  FAC 

SNS  (2)  = ZTS  FAC 

RETURN  FAC 

END  FAC 
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Appendix  B 

Block  Diagram  and  Flow  Chart  for  WIMP 
Ray-Tracing  Program  TRISL 


Explanation: 

HI.  TRISL  BLOCK  DIAGRAM 

This  presentation  shows  the  general  communication  network  among  the  various 
functional  processes  of  TRISL,  and  illustrates  the  iterative  complications  inherent 


in  the  program. 

The  function  of  the  various  blocks  is  generally  described  as  follows 

BL.OCK  I: 

Initialization  Procedures. 

BLOCK  II: 

E-layer  refraction  upwards. 

BLOCK  III: 

Fj -layer  refraction  upwards. 

BLOCK  TV: 

Reflection  layer  initialization. 

BLOCK  V: 

Tilt  calculation  and  consistency  checks. 

BLOCK  VI: 

Fj -layer  refraction  downwards. 

BLOCK  VII: 

E-layer  refraction  downwards. 

BLOCK  VIII: 

Error  checks. 

BLOCK  IX: 

Iterate  reflection  and  tilt. 

BLOCK  X: 

T raffic  control. 

BLOCK  XI: 

Normal  Terminal  calculations. 

BLOCK  XII: 

ERROR  exits. 
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112.  DETAILED  FLOWCHARTS 


Detailed  flow  charts  are  presented  for  each  block.  Most  line  members  of 
TRISL  are  accounted  for  except  when  so  doing  unduly  complicates  comprehension 
of  the  process. 


H3.  FLOW  CHART  SYMBOLS  AND  CONVENTIONS 


(a)  Process  entry/exit  point. 

(b)  Processing  statements. 

(c)  Decision  trees,  two  or  three  way  branches. 

(d)  External  procedures. 

(e)  Numbers  in  the  form  Snnn  refer  to  statement  numbers  in  the  program. 

(f)  Numbers  in  the  form  nnn  refer  to  line  numbers  of  the  program  listing. 

(g)  Entry  points  to  a block  are  labelled  generally  with  block  designation  state- 
ment number  and  line  number  specifying  the  .source. 

(h)  Exit  points  from  a block  are  generally  labelled  by  block  designation,  state- 
ment number,  and  line  number  specifying  the  destination. 

(i)  Error  exits  are  labelled  by  an  error  number  and  the  error  condition 
detected. 

(j)  Exits  from  the  bottom  generally  enter  the  next  block  at  the  top. 

(k)  External  call  names  and  calling  line  numbers  are  contained  within  the 
external  procedure  symbols. 


IU.  PROGRAM  LISTING 

Listings  are  included  for  TRISL  together  with  its  externals. 
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TRISL  Block  Diagram 
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TRISL:  Block  1,  Argument  List 


enient 

1/D 

Mnemonic 

Type 

1 

I 

MAI 

Double 

2 

I 

NLA  I 

Double 

3 

1 

W LAI 

Double 

4 

I 

NLAZ 

Double 

5 

I 

WLAZ 

Double 

6 

I 

MAZ 

Double 

7 

I 

FREQ 

Double 

8 

I 

MODE 

Integer 

9 

I 

RTARG 

Double 

10 

O 

RS  2 

Double 

11 

O 

nltarg 

Double 

12 

o 

WLTARG 

Double 

13 

o 

PPTOT 

Double 

14 

o 

GPTOT 

Double 

15 

o 

AZF 

Double 

16 

o 

ELF 

Double 

17 

o 

SRANGE 

Double 

18 

o 

ELSR 

Doub  le 

19 

o 

DIST 

Double 

20 

o 

BEAR 

Double 

21 

o 

IGOBAK 

Integer 

Definition 

Geocentric  Distance,  Initial  Point 

North  Latitude,  Initial  Point 

West  Longitude,  Initial  Point 

Takeoff  Bearing  (all  calls)  or 
Target  North  Latitude 

Takeoff  Elevation  (all  calls)  or 
Target  West  Longitude 

Less  than  6000  km  (all  calls)  or 
Target  Geocentric  Distance 

Frequency 

| MODE  I = Number  of  hops 

= 6370  km,  (all  calls) 

Geocentric  Distance,  End  Point 

North  Latitude,  End  Point 

West  Longitude,  End  Point 

Phase  Path  Length 

Group  Path  Length 

Final  to  Initial  Azimuth 

End  Point  Elevation 

Final-Initial  Straight  Line  Distance 

Takeoff  Elevation 

Final -Initial  Great  Circle  Ground 
Range 

Initial  to  Final  Azimuth 
Error  Flag  8 = 0 No  Error 

- 1 Fatal  Error 
= 2 Penetrate 
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BLOCK 

xn 


ALLOCATE  STORAGE 
SPECIFICATIONS 


INITIALIZE  VARIABLES 


1 

— 

n 

77- 

- 85 

87-102 

INPUT  PARAMETERS 
CORRESPOND  TO  ELEVATION 
AND  AZIMUTH.  ONLY  OPTION 
USED  IN  ALL  DRIVING 
PROGRAMS 


S90 


INITIALIZE  TAKEOFF 
VECTORS 


FROM 

BLOCK 

X 


NOTES 

NUMBERS  IN  PROCESSING 
SYMBOLS  REFER  TO 
PROGRAM  LINE  NUMBERS. 

Snnn  REFERS  TO  FORTRAN 
STATEMENT  LABEL  nnn. 


Block  I,  Initialization  Procedures 


61 


FROM  BLOCK  I 


BLOCK  /S59I 
XI  ! 1049 


125-157 


160-173 


E- LAYER 
REFLECTION 

BLOCK 

nz  Vi96y 


1 175- 

Cl 

SI  30  ] 

tz 

| 182- 

-192  J 

\BLOCK 

487  J V. 


f*  BLOCK 
1 906  J IX 


xn 

^rr 

BLOCK 

/s200\- 

REFRACTION 

m 

\267j^ 

SOLUTION 

BLOCK 

/S2I5\ 

CONVERGE 

m 

I*™ ; 

FAILURE 

211-2271 

<^228> 


239-252 


254-261 


BLOCK  f ERR 

xn  \S TF 


TOO  MANY 
ITERATIONS 


Block  II,  E- Layer  Refraction  Upwards 
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Block  V,  Tilt  Calculation  and  Consistency  Checks 


65 


66 


741  -748 


750-752 
I STOP =8 


749  ^ 

<3^ 


753-770 


789-791 

IST0P=9 


773-775 


S4305 


776-781 


783-785 

IST0P=I5 


1806-8081 

IST0P=I5 


S43 


803-804 


809-824 

1 STOP  = 1 

A 

I STOP  = I 


ISTOP=l 


871-891 


INCONSISTENCIES  DETECTED 
TO  BE  RESOLVED 


/S475> 

\ BLOCK 

/S50i ' 

\ BLOCK 

/S5iO\ 

^836, 

/ IX 

\ 893  > 

‘ IX 

V923 ) 

NOTES 

BLOCK  X EXIT  IS  DESIRED  CONVERGENT  SOLUTION  TO 
TILTED  LAYER  ITERATION 

THE  ISTOP=  I EXITS  TO  S475  ARE  NORMAL. THE  OTHER 
ERROR  CONDITIONS  MAY  BE  CORRECTED  IN  SUCCEEDING 
ITERATIONS 

Block  VIII,  Error  Checks 


FROM 

BLOCKS  n m UL 


Block  XI,  Normal  Terminal  Calculations 


OUTINE  TRISL 


7 4/74  OPT  = l 


FTN  4.5*414 


04/ 27/ 


SUBROUTINE  T RISK  HAi,  NL  A1 , WLA1 , NLA2  , HLA2  , HA2  , FREQ, MODE  ,R  TARG  ,RS  2 , TR2 
*NL TARG, WL TAR G,PPTOT ,GPTOT  , AZF , ELF, S RANGE ,ELSR,DIST,BEAR, IGOBAK)  TR2 
IMPLICIT  00U8LE  PRECISION  (A-H.O-Z)  TR2 

LOGICAL  HITS, FLAG, SEARCH, A SC,  FLAG2,FLAG3, NOEL  AY, INCHOP, LP<  6)  , NIP  TR2 

LOGICAL  GPONLY  TR2 

DOUBLE  PRECISION  MESSPR  TR2 

DOUBLE  PRECISION  NL A1 ,HA 1, NL A2, M A2, NLT ARG, HAG , MA1B1  , HBi, NLC1 , NL4 , TR2 
1MRX1SQ, NLC2,  MR3  , NL8  , NLC3 ,MNi,HA2C2 ,MC2,MA3C3 ,MC3 ,NLD ,HRH 1 , HRH2  , TR2 

2HRY 1SQ ,MC4S0 ,NLE , MRH3 ,MRH4 , MC5 SO , NLH2 ,KP, L, Ni , H,  KPPP,LOS,  TR2 

3H  AG  VEC,  NORTH  ( 3)  , ME,  MN,  NORTH  i(  3)  ,E  ASTI  (3  > , HR2  , MCI  PSQ  , H3  00  0 TR2 

REAL  0UT2  TR2 

DIMENSION  Sl(3) ,Ai(3>  ,Bi(3> , A 11 ( 3) , E ASTC3) ,MESSPR(  15) , TR2 

* XI(3),KP(3)  , L (3 ,3) *P1P(3)»PC2(3),P8(3)»R3(3),PC3(3)»P10(3),C  TR2 

*2(3),PRl(3),Ni(3)*P0(3)»A2(3),M(3,3) , C3(3),A3(3)  , PH  1(3), KPPP(3) ,P(  TR2 
*3,3)  ,PD(3I  ,PE(3)  ,PQ(3)  ,PW2(  3)  ,W2(3)  ,PTARG(3)  , RN1(3),R2(3  TR2 

*) ,PC1 (3) , P4(3)  TR2 

COMMON/ FLIHSY/0IST1,I CNPRB  /GPFLA G/GPONLY ,IHRITE  T R2 

COMMON/CPREY/P0,  P,RTi,C52PPP,C53PPP, MRH1,CBETM1,HRM4,RY,  TR2 

2 0PP1P,MC5SQ,RM1, YY2PPP, ZY2PPP,H1,MC4SQ,MRY1SQ,YY1PPP,  TR2 

2 ZY1PPP ,C42PPP  »C43PPP ,YS 3PPP  , ZS  3PPP, C32PPP.RTM1,  ONDR, DNDNL , TR2 

3 0 NO ML, 0, GNU,  PP2,  PPEF 1D,CB ETH2, H E, PPF 2,RV,HF 1,PP2A ,DPP 2P , TR2 

4 MBl,SBETi,GPF2,GPEFiO,OGP2P,DGPiP,  T R2 

5 FI, PREV(56) , HITS, INCHOP, FLAG,  FLAG2, FLAG3,N0ELAY,LP  TR2 

COMMON/ VEC TRS/LOC (3, 8), VEC (3, 6) ,IREFL  TR2 

CONMON/NIPRIP/NIP  TR2 

•/RIIPAR/H300  0 , FCF2  , FCFi , FCE , HBE ,HAE, HME, HAF1 , HHF 1,HAF2»HMF2,DUM(13  TR2 
*) ,10(3) /RTCOM/COSBO,NUSE  TR2 

COMMON/REMWEN/INUSE,NNUSE( 20) ,FCREF(20>  TR2 

DATA  B, EPS1.EPS 2, EPS3 ,EPS4 , EPS5/6670 .DO ,5*1 .D  0/  TR2 

DATA  PI, RTO, OT R, RO/3. 1415 92 653589793  DO, 57. 29 577951308 23  2D  0,  TR2 

*.0174532925199433000, 6370. DO/  TR2 

DATA  LIMEA,LIMF1A,LIMF10, LIMED, LIMF2/5*15/  TR2 

0  AT  A MESSPR/  7HF2  ITER,8HBIG  TILT,8HBIG  TILT  ,9HF1-D  ITER,  TR2 

1 9HF1-D  REFL.7HF1  MISS,8HE-D  ITER, 8HE-D  REFL, 6HE  MISS,  TR2 

2 8HI-I  REFL, 9HE-F 1 OUCT,8HE-A  ITER, 10HF1-F2  DUCT,  TR2 

3 9HF1-A  ITER,9HMISS  TARG/  TR2 

MLON  ( ARG1 , ARG2)  = DSIGN(PI,ARG1) +PI-DATAN2(ARG1, ARG2)  TR2 

0 ACOS ( ARG 1) =0  AT  AN2 (OS  QRT  C1.D0-ARG1**2),ARG1)  TR2 

DASIN(ARG1)=OATAN2(ARG1,OSQRT(1.DO-ARG1**2))  TR2 

INCHOP= • TRUE . TR2 

HITS=.TRUE.  TR2 

ASC=MODE • EQ.O  TR2 

SEARCH=ASC  TR2 

PPO=O.DO  TR2 

PPEF1A=  0.00  TR2 

OPP1=0. DO  TR2 

PP2  A=0.  DO  TR2 

PPF  2=0.00  TR2 

PP2=0 ,D0  TR2 

PPE  FI  0=  0 . DO  TR2 

0 PP  2P=0 • DO  TR2 

PP0P=0.O0  TR2 

PPT OT  =0 .0  0 TR2 

GPT  OT=0 .00  TR2 

GPE  Fi A=0 • DO  TR2 

DGP  1=0.  DO  TR2 

DGP  2=0.00  TR2 
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74  m OPT  = 1 


FTN  4.5+414 


04/27/ 


GPEF1A=  0.00  TR2 

GPF2=0.00  TR2 

GPEF1D=0.D0  TR2 

DGP2P=D.D0  TR2 

B GP  .0  0 TR2 

1CNPRB=0  TR2 

INiiSE=0  TR2 

IG08AK=  0 TR2 

IHOPS=  0 TR2 

RS2=RTARG  TR2 

CN=OCOS  (NLA1+OTR)  TR2 

SN*OSIN (NLA1 +OTR)  TR2 

CW=DCOS (WLA1*0TR)  TR2 

SH=DSIN (WLA1+OTR)  TR2 

Alt  1) =CN*CW  TR2 

A1C2)=-CN*SW  TR2 

A1(3)=SN  T R2 

IFCMA2.LT. 6000.00)  GO  TO  50  TR2 

CN=0C0S(NLA2*0TR)  TR2 

DO  45  1=1,3  TR2 

A11(I)=MA1*A1(I)  TR2 

45  A1  ( I)  =A11  ( I)  TR2 

Bit  ll3MA2+DCOS(WLA2*OTR)  +CN-A1C1)  TR2 

B1C2)=-MA2+OSIN(HLA2+OTR)*CN-A1 C2)  TR2 

81(3)  = HA2+DSIN(  NLA2  *OTR)  -A1C  3)  TR2 

MB1=H AG (Bll  TR2 

OTA1B1=OOT(A1,B1)  TR2 

GO  TO  60  TR2 

C WHEN  HA2  IS  INPUTTEO  AS  A NUMBER  LESS  THAN  6000,  NLA2  IS  ASSUMED  TR2 

C TO  BE  THE  BEARING, AND  HLA2  IS  ASSUMED  TO  BE  THE  TAKE-OFF  ANGLE  TR2 

50  EAST 1 (1)=SW  TR2 

EAST1 (2)=CW  TR2 

E ASTI (3 ) = 0 • DO  TR2 

CALL  CROSS (A1,EAST1,N0RTH1)  TR2 

CAZ=OCOS(NLA 2+OTR)  TR2 

SAZ=0SIN(NLA2»0TR)  TR2 

CEL=OCOS( HLA2+OTR)  TR2 

SEL=DSIN(HLA2*0TR)  TR2 

DO  55  1=1,3  TR2 

B1 (I) =C  EL* (CA  Z’NORT  Hi (I ) +SAZ*EAST1 Cl ) )+SEL* A1 (I)  TR2 

All (I)= MA1+A1 ( I)  TR2 

55  AKDaAllCI)  TR2 

HB1*1.D0  TR2 

OTA  IB 1* SEL*MA1  TR2 

60  IF (DABS (MA1-R0) .LE.l.D-5)  GO  TO  90  TR2 

IF(OASIN<RO/MA1)-OACOS<-OTA1B1/(MA1*HB1>) > 86,88,80  TR2 

80  T 1=  (-DT  A1 B1-OS3RT CDTAiBi+*2-DOT  (B1 , Bl)  * (DOT  (A1,A1)-RQ**2>) ) TR2 

*/OOT(Bl,Bl)  TR2 

DO  85  J=1 ,3  T R2 

85  A1(J)=B1(J)*T1+A1(J)  TR2 

A=-2.DO*DOT ( A 1 , B 1 ) /DOT ( A1 , A1 ) TR2 

00  87  J*l,3  TR2 

87  81(J>=B1(J)+A*A1(J)  TR2 

PPT0T=T1*MB1  T R2 

GPT  OT=PPT  OT  TR2 

GO  TO  90  TR2 

88  IF(OTA1B1.GE.O.OO)GO  TO  90  TR2 
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PPTOT=-2.O0*DTAlBl/H81  TR2 

GPT  OT =PPTOT  TR2 

90  IF C INCHOP)  GO  TO  92  TR2 

IONREF= IONREF+1  TR2 

IF (IONREF .LT .40)  GO  TO  93  TR2 

1ST  OP=10  TR2 

GO  TO  7310  TR 2 

92  I ONREF=  0 TR2 

IHOPS=IHOPS+l  TR2 

93  PREV( 13) =0 .00  T R2 

C START  OF  E LAYER  REFRACTION  TR2 

NIP=. FALSE.  TR2 

INT=0  TR2 

INUSE= INUSE+1  TR2 

FLAG=. FALSE.  TR2 

FLAG2=. FALSE.  TR2 

CALL  CROSS( Ai,Bl,XI)  TR2 

CALL  CROSS ( XI i B li KP ) TR2 

MAI Bl =MAG ( XI)  TR2 

MB1=MAG (Bl)  TR2 

A 1MAG=M AG ( Al)  TR2 

C EON (2)  TR2 

DO  100  J= 1 i 3 TR2 

L(1,J)=XI(J) /MA1B1  TR2 

L (2,J)=B1( J)/M31  TR2 

100  L<3,  J)  = KP(J)/(NA181*MB1>  TR2 

Y5  P=0  • 00  TR2 

Z5P=0.00  TR2 

C EON (9 ) TR2 

PI  P (1  )=  0.  DO  TR2 

DO  110  1=2,3  TR2 

110  P1P(I)  = L(I,1)*A1(1)  +L(I,2)  *A1(2) +L(I,3) *A1(3)  TR2 

YC1P=DSQRT(B**2-P1P(3)**2>  TR2 

C EQN (12)  TR2 

DO  120  J=  1 , 3 TR2 

120  PC1(J)=L(2,J)*YC1P+L(3,J)*P1P(3>  TR2 

C EQN (14, 15)  TR  2 

NLC1=0A  SIN (PC  1 ( 3 ) / B ) TR2 

MLC1=WL0N(PC1  (2)  , PC1(  1)  ) TR2 

C EQN  ( 17)  TR2 

BETACl=DASIN(OOT (PCI,  Bl)  / ( MB 1+B) ) TR2 

8CBC1=B*0C0S(BETAC1)  TR2 

125  CALL  RIIP(R0,NLC1, HLC1 , PL A SO)  T R2 

I F< SEARCH. AND.  ASC.  ANO.  RT ARG . GT.BCBC1. AND. RTA RG.LE.RO +H BE) GOTO  591  TR2 
NIP=  • TRUE • TR  2 

NOELAY=. FALSE.  TR2 

IF  (FCE.GT .O.DO) GO  TO  126  TR2 

FCE=l.D-2  TR2 

NOELAY  = .TRU£.  TR2 

126  RH=RO+HBE+HAE  TR2 

H=H AE  TR2 

F=FREQ/FCE  TR2 

Y3P=0SQRT(RM**2-P1P(3)**Z>  T R2 

Y8EP=0SQRT ( (RO  + HBE)  *»2-PlP(3) **2)  TR2 

DO  1265  J=l, 3 TR2 

1265  PC1(J)=L(2,J)*Y3P+L(3,J)*P1P(3)  TR2 

NLClr  DASIN(PC1  ( 3)  /RM)  TR2 
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WLC1=WL0N(PC1  (2)  , PC  1(1)1 
CALL  RIIP(R0,NLC1,WLC1  ,PLASD) 

127  RM=H9E*HAE+R0 
H=H  AE 

NOELAY=. FALSE. 

IF(FCE.GT.O.OO)GO  TO  129 

FCE=l.D-2 

NOELA Y= .TRUE • 

12  9 F=FREQ/FCE 

130  IF  < .NOT .FL  AG) GO  TO  1305 
CBET1=  DA8S(  Z 10PP)  /RH 
SBET1  = DSQRT(1.00-C BET  1**2) 

Y 3P=RH*  SBET1 
YBEP=(R0+HBE)*SBET1 

GO  TO  131 

1305  Y3P=DS0RT(RM**2-P1P (3)**2) 

YBEP=  DSQRT  ( ( RO  *HB  E)  **2-PlP(3>  **2) 

C EQN(lfl) 

C BETi  = BCBC 1/RM 

S8ET1=DSQRT(  1. OO-C BET 1*C  BET1) 

131  IF(F*SBET1.GT.1.00)  GO  TO  135 
IF( FLAG)  GO  TO  134 

132  IF (FL AG2)  GO  TO  134 
1325  IREFL=1 

RM=R0*H9E*HAE 
HR  3 = A 1M  AG 

SIN02=D0T(Al,Bl)  /(MB1*A1HAG) 

COSO2  = OSQRT(l.D0-SI  N02**2 ) 

Y9P=P1P ( 2) 

Z 9P=P IP  (3) 

C2  2P=L( 2,1) *01(1)  *L (2,2)*B1(2) *L (2, 3) *B 1 ( 3) 

C23P=L( 3, 11*81(1)  *L <3,2)*B1(2)*L(3,3)*B1(3> 

Y10P=COSO2*( Y9P*C0S02*Z9P*SIN02> 
Z10P=COSO2*(Z9P*COSO2-Y9P*SINO2) 

GO  TO  3 096 
134  FLAG=. FALSE. 

1ST  OP=l 1 
GO  TO  7310 

135  IF ( .NOT .NOELAY)  GO  TO  137 
Y4P= Y3P 
OGP1=O.QO 
Z 4P=P IP  (3) 

GO  TO  139 

13  7 OGPl=-H/SBETl  + H*F/2.O0*OLOG(  (F*  SBET1 +1  .00)/ (F*SBET  1-1.D0)  ) 

GAH1=CBET1*0GP1/RM 

Y4P=Y  3P*OCOS (GAM1)-P1P(3)*DSIN(GAH1) 

Z4P=P1P (3) *OCOS (G AMI)  *Y3P*DSIN(GAM1) 

C DIST  BETWEEN  SUCCESSIVE  R<S 

13  9 OISTB=DSQRT(  (Y4P-Y5P)  **2*(  Z4P-Z5P)  **2> 

IF ( FLAG ) GO  TO  141 
00  140  J=l,3 

140  P4(J)=L(2,J)  *Y4P+L(3,J)*Z4P 
GO  TO  1418 

141  00  1413  J=1 , 3 

1413  P4(J)=M  (2,J)*Y4P*H(3»J)*Z4P+P0(J) 

1418  IF (01 ST8.L  E. EPS1)  GO  TO  200 

IF( INT. LE.l. OR.OISTB.LT.  OISTA ) GO  TO  145 
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0UT2=0ISTA  TR2 

1F(  IWRITE.EQ.i) WRITE<6,142>  OUT2.IHOPS  TR2 

142  FORMAT ( 24H  E-LAYER  CONV.  PROBLEMS, , F6. 1, 17H  KH  USED,  HOP  NO. ,13)  T R2 

H=HP  TR2 

Y3P=T3PP  TR2 

RM=  RMP  TR2 

DGP 1=GPP  TR2 

F=FP  TR2 

GO  TO  215  T R2 

145  NL4=0 AS INCP4(3) /MAG(P4) ) TR2 

WL4=WL0NtP4(2),P4<l))  TR2 

DISTA=0 ISTB  T R2 

GPP=DGP 1 TR2 

V5P=Y4P  TR2 

75P=Z4P  TR2 

HP=H  TR2 

FP=F  TR2 

Y3PP=Y3P  TR2 

RMP=RM  TR2 

SIN01=S8ET 1 TR2 

C0S01=C  BET1  TP.2 

DO  148  J=  1 , 3 TR2 

148  R2(J)=P4(J)  TR2 

CALL  RIIP(RO,NL4,HL4,PLASO)  TR2 

NOELAY=. FALSE.  TR2 

IF(FCE.GT.O.DQ)  GO  TO  149  TR2 

FCE-1 .0-2  T R2 

NOEL AY= . T RUE • TR2 

149  RM=  HB  E+HA  E*R  0 TR2 

H=HAE  T R2 

F=FREO/FCE  T R2 

I NT  = I NT  *1  TR2 

IF(INT.LE.LIHEA)  GO  TO  130  TR2 

IST0P=12  TR2 

GO  TO  7310  TR2 

C ENO  OF  E LAYER  REFRACTION  T R2 

C START  OF  FI  LAYER  REFRACTION  TR2 

200  Y5P=Y4P  T R2 

Z5P=74P  TR2 

SIN01=SBET 1 TR2 

C0S01=C  BET  1 TR2 

00  210  J = 1 , 3 T R2 

210  R2(J)=P4(J)  TR2 

215  MR  2=RM  TR2 

RR1-HR2-H  TR2 

SINOl  P=  SINOl  TR2 

COS01P=COS01  T R2 

Y 6P=COS 01* <Y5P*C0S01+Z5P* SINOl)  TR2 

Z6P=COSOl*<Z5P*COS01-Y5P*SINOl>  TR2 

IF( .NOT. FLAG)  GO  TO  2155  TR2 

00  2151  J = 1 , 3 T R2 

2151  P4(J)=M(2,J)  *Y3P*H(3,  J)*Z10PP  TR2 

RR1=-H/MAG(P4)  TR2 

00  2152  J=1 , 3 TR2 

2152  P4(J)=R2(J)*RR1*P4(J)  TR2 

RR1=M AG (P4)  TR2 

MR2=MAG(R2)  TR2 
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COSOiP=BCBCl/MR2  TR2 

SIN01P=0SQRT( l.DO-COS01P**2>  T R2 

00  2153  J=  1 , 3 TR2 

P4(J)=M(2,  J)  + Y6P+M(3,  J)*Z6P+P0(J)  TR 2 

2153  P8( J) =R2( J) -P4( J)  TR2 

CALL  CROSS (R2, P8,  XI)  TR2 

CALL  CR0SS(XI,P8, KP)  T R2 

CNST=MAG(XI)  TR2 

FL0G=MAG(P8)  TR2 

00  2154  J=  1 , 3 TR2 

L (1  » J)  = XI  (J)  / CNST  TR2 

L (2, J)  = P8 ( J) /FLOG  TR2 

2154  L(3,J)  = KP(J)/ (CNST*FLOG)  TR2 

Y5P=L(2,1)*R2(1)+L(2, 2) *R2 (2) +L (2 ,3) *R2 (3 ) TR2 

Z5P=L  (3 ,1>*R2  (1)  +L  (3,  2)  *R2(2)  +L(  3,  3)  *R2(3)  TR2 

Y6P=L(2 ,1)  *P4(1) +L (2,2) *P4(2>  +L (2, 3)*P4(3)  T R2 

Z6P=L (3,1) *P4 (1) +L(3,2) *P4 (2) +L(3,  3)  *P4(  3)  T R2 

FLAG=. FALSE.  TR2 

2155  INT  =0  TR2 

FLAG3=. FALSE.  TR2 

HE  A=H  TR2 

Hc  = H TR2 

'■’PEF 1 A=  0.  0 0 TR2 

IF (.NOT. NOEL A Y)PPEF1A=H* (.5O0  + SINO 1-1. 0 0/S INOl+F/4. DO*  (1.00-  F** (-2  TR2 

*) +C0S01**2)*0L3G( (F’SINO  1 + 1.00) / (F*SIN01-1. DO ) ) ) TR2 

3PP1=  PPEF1 A TR2 

>PO=RR1*OSQRT(1.00-  (C  OS  01  P*  HR2 /RR1 ) * *2  ) TR2 

1 AU1EA=MR2*SIN01P-PPO  TR2 

IF(. NOT. HITS. OR. COS01P*MR2/A1MAG.GE. 1.00)  GO  TO  217  TR2 

PPO=PPO-A1MAG*OSORT  (1.0  0-(COS01P*MR2/A1MAG) **2)  TR2 

217  Cl  2 P=  Y5  P-  Y5  P T R2 

C13P=Z5P-Z6P  TR2 

MC1PSQ=C12P**2+C13P**2  TR2 

HRXlSa=BCBCl**2  TR2 

218  RM=R0+HBE+HAE+HAF1  TR2 

TC2=0SQRT((RM**2-MRX1SQ) /MC1PSQ)  TR2 

IF( SEARCH.  ANO.RTARG  ,GT. OSQRT (HRX1SQ)  .ANO.RTARG.L  E.RR1)  GO  TO  591  TR2 

YC2P=C12P*TC2+Y6P  TR2 

ZC2P=C13P*TC2+Z6P  TR2 

DO  220  J= 1 ,3  TR2 

220  PC2  (J)=L(2,  J)*YC2P+L(3,  J)  *ZC2P  T R2 

NLC2=  OAS IN (PC2 (3) /RM)  T R2 

HLC2  = HLON  (PC2  (2) , PC2(  1)  ) TR2 

CALL  PUT  (RO,NLC2,WLC2,PLASO)  TR2 

225  RM=R0+HBE-'HAE  + HAF1  TR2 

H=H AF1  TR2 

F=FREQ/FCF1  TR2 

Y9P=  0.0  0 TR2 

Z9P=0.D0  TR2 

230  IFf  .NOT. FLAG)  GO  TO  2305  TR2 

C BF. T 4 = OABS  (Z10PP)/RM  T R2 

S8ET4=0SQRT(  1 . DO -C BET4* *2 ) TR2 

Y 7P  = RM*  S8ET4  TR2 

GO  TO  231  TR2 

2305  T0EL3=DS0RT( (RM**2-HPX1SQ) /MC1PSQ ) TR2 

Y7P=C12  P*  T OE  L3  + Y6  P T R2 

Z7P=C13P*T0EL3+Z6P  TR  2 
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NR3  =RM  TR2 

CBET4=BCBC1/MR3  T R2 

S BET  4=DSQ  RT(  1. DO- CBET4+CBET4  ) T R2 

231  IF(F+SBET4.GT.1.00»  GO  TO  235  TR2 

IF (FLAG)  GO  TO  234  TR2 

232  IF<FLAG3>  GO  TO  234  TR2 

2325  IRE  FL=2  T R2 

RM=R0+HRE+HAE+HAF1  TR2 

Y9P=Y5P  TR2 

Z9P=Z5P  TR2 

MR3=MR2  TR2 

COS02=COS01  TR2 

SIN02  = S INOl  TR2 

Y10P=COSO2*< Y9P*C0S02  + Z9P*SIN02  > TR2 

710P=COSO2MZ9P*COS02-Y9P*SINO2)  TR2 

GO  TO  3095  TR2 

234  FLAG=. FALSE.  TR2 

I ST  OP=i 3 TR2 

GO  TO  7310  TR2 

23  5 OGP2=H*  (F/2. 0 0*0 LOG  ((F+SBET4+1.  OO  » /(F*SBET4-1  .00  ) ) -1.D0/SBET4)  TR2 

GAM2=CBET4*0GP2/RH  TR2 

Y8P= Y7P*0C0S(GAM2) -Z7P*0S I N ( G AM  2)  TR2 

ZBP=Z7P*DC0S  (GAM2)  +Y7  P*  OSI N ( GA  M 2)  T R2 

C 0 1ST  BETWEEN  SUCCESSIVE  R<S  T R2 

DIST8  = DSQRT( ( Y 8 P- Y9P) **2 + (Z 8P-Z 9P)  ** 2 ) TR2 

IF(FLAG)  GO  TO  241  TR2 

DO  240  J=i,3  TR2 

240  P8(J)  = L (2, J) *Y3P+L(3, J)  *Z8P  T R2 

GO  TO  2425  TR2 

241  no  242  J=l,3  TR2 

242  P8(J)=M(2,  J)  *YBP  + M<3,  J)  *Z8P  + P0(J>  TR2 

2425  IF(OISTB.LE.EPS2)  GO  TO  300  TR2 

IF(INT.LE.1.0R.0ISTB.LT.DISTA)G0  TO  245  TR2 

OUT2=  01  ST  A TR2 

IF(IWRITE.EQ.l)  WRITE!  6,  24  3)  OUT2rIHOPS  TR2 

243  F ORMA  T ( 25  H PI -LAYER  C ONV . PRO BL EM S , ,F 6 . 1 , 1 7H  KM  USED,  HOP  NO., 13)  T R2 

H=  HP  TR2 

F=FP  TR2 

Y7P=Y7PP  TR2 

RM=RMP  TR2 

OGP2=GPP  TR2 

GO  TO  307  TR2 

245  NL8=  DASIN(P8(3)/MAG(P8>)  TR2 

WL8=WLON(P8(2),P8(l))  TR2 

Y9P=Y8P  TR2 

Z9P=  78P  TR2 

HP  = H TR2 

FP=  F TR2 

Y 7PP=  Y7  P TR2 

RMP=RM  TR2 

GPP  = DGP  2 TR2 

COSO  2=  f RE T 4 TR2 

S IN02=S  BET  4 TR2 

DO  247  J= 1 , 3 TR2 

247  R3  ! J ) =P  8 ( J ) TR2 

PISTA=n ISTB  TR2 

CALL  RIIP(R0,NL8,Wl8,PL4Sn>  TR2 
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RM=  R0+H8E+HAE+HAF1  TR2 

H=H  AF  1 TR2 

F=FREQ/FCF 1 T R2 

I NT  = I NT +1  T R2 

IFIINT. LE. LIMF1 A)  GO  TO  230  T R2 

IST0P  = 14  T R2 

GO  TO  7310  TR2 

C ENO  OF  FI  LAYER  REFRACTION  T R2 

300  Y9P=Y8P  TR2 

Z9P=Z0P  T R2 

SIN02=S8ET4  TR2 

C0S02=C  8ET4  T R2 

00  305  J=l,3  TR2 

305  R3(J»=P8<J)  TR2 

307  MR3=RM  TR2 

IREFL=3  TR2 

RM=  RO  + HHF2  T R2 

RR2=MR3-H  TR2 

Y10P=COS02* <Y9P*C0S02*Z9P*SIN02)  TR2 

Z10  P=C0S02  *(  Z9P*C0S02-Y9P*SIN0  2)  T R2 

C0SO2P=COSO2  TR2 

SIN02P=SIN02  TR2 

IF(. NOT. FLAG)  GO  TO  3075  TR2 

DO  3071  J=l,3  T R2 

3071  P8(J)=M(2, JJ *Y7P*M(3,J)*Z10P»PO(J)  T R2 

RR2=-H/MAG(P0)  T R2 

DO  3072  J=l, 3 TR2 

3072  P8 ( J ) =R3 (J)*RR2*P8 (J)  T R2 

RR2=MAG  ( P8)  TR2 

HR3=M  AG (R  3)  TR2 

C0S02P=  0CBC1 /MR3  TR2 

SIN02P=0SQRT (i.D0-COS02P**2)  TR2 

DO  3073  J=l, 3 TR2 

P8(J)=M«2,J)  *Y10P*M(3,J)*Z10P*P0(J)  TR2 

3073  P4(J)=R3(J)-P8(J)  TR2 

CALL  CROSS!  R3  , P4  , XI)  T R2 

CALL  CROSS(XI,P4,KP)  T R2 

CNST=MAG ( XI)  TR2 

FL0G=NAG(P4)  TR2 

00  3074  J=l, 3 TR2 

Ltl , J)  = XI(  J) /CNST  TR2 

L(2,J)  = P4  (J)/FLOG  TR2 

3074  L(3,J)=KP(J) /<CNST*FL0G>  TR2 

Y9P=L(2,1)*R3(1)*L<2,  2)  *R3(  2)  *L  ( 2,  3)  *R3(3)  TR2 

Z9P=L(3,1>  *R3(1) +L  (3,2)*R3(2> +L ( 3 , 3)  * R3 < 3 ) T R2 

Y10P=L<2,1  )*P8  (1 ) *L  (2, 2 ) *P8  (2)  +L  C 2,  3)  *P  8<  3)  T R2 

Z10P=L(3,1)*P8(1)+L(3,2)*P8(2)*L(3,3)*P8(3)  T R2 

FLAG=. FALSE.  TR2 

3075  DO  308  J=l,3  TR2 

308  PREV(J)=0.00  TR2 

IF(. NOT. SEARCH)  GO  TO  309  T R2 

IF(RT  ARG.  GT . RR1. AND. RTARG.  LE.  HR2 ) GO  TO  593  T R2 

IF (IREFL • EQ.  1)  GO  TO  309  TR2 

IF1MR2.LE.RR2.AN0. RTARG. LE.RR2. AND. RTARG. GT.MR2)  GO  TO  593  TR2 

IF(RR1.GT.MR3. AND. RTARG. LE.RR1. AND. RTARG. GT.HR3)  GO  TO  593  TR2 

C TILTED  F2  LAYER  TR2 

309  HF 1 A=H  TR2. 
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HE  1=H  TR2 

TAU=MR3*SIN02P-RR2+0S0RT < 1.  D 0- (M R3/RR2+ COS  02  P ) * + 2 ) + T A U1E  A TR2 

PPE  FI  A=PPEFl  A +H*  ( .5DO*SIN02+F/4.DO*(  1.0  0-F**(  - 2)  + C0S0  2»*2)  *DLOG  ( (F  T R2 
+ + SINO2+1.D0) /(F+SINO2-i.OQ))-i.D0/SINO2>  +TAU  TR2 

GPEF1A=0GP1+0GP2+TAU  TR2 

3095  C22P=Y9P-Y10P  TR2 

C 23P=  Z9P- 710P  TR2 

30  96  TC3=(RM*+2-Y10P++2-710P**2) / (C2 2P **2 +C 23P ** 2 ) T R2 

IF(TC3.  LT.  0 . OO ) GO  TO  (134 , 134, 234) , IREFL  T R2 

TC3=OSOPT (TC3)  TR2 

YC3P=C22P*TC3+Y10P  TRP 

ZC3P=C23P*TC3+Z10P  TR2 

00  310  « 1,3  r R2 

PC  3 (J)  = L (2, J)*YC3P  + L(3, J) +ZC3P  TR2 

P10(  J)=L(  2,  J)  *Y10P+L(  3,J)  *Z10P  TR2 

310  C2 (J)=L (2, J)*C22P+L (3  ,J)*C 23P  TR2 

BET  A22=0 . DO  TP2 

NLC3=0ASIN(PC3(3) /RM)  TR2 

WLC3  = WL  ON (PC  3 ( 2) , PC3 ( 1) ) TR2 

N0FLAG=  0 TR2 

CALL  HTGROF(  NOFLAG,  FREQ,  R TD’NLC  3 , RTD  *WL  C3 , RTD  *DAC0S  ( BC  BC  1/RM ) ,0.  T R2 

*0  0 , RM , P Ml ,RT 1 ,RTM1  , ONOR, ONONL , DNDW L , PNTFLG)  TR2 

NOEL AG= 1 TR2 

GO  TO  (311,312,313) , IREFL  TR2 

311  F1  = FREQ/FCE  T R2 

H1=HAE  TR2 

I F ( NUSE • EO • 1 ) GO  TO  314  TR2 

FL A G2=. TR  UE.  TR2 

INT  = 0 TR2 

GO  TO  127  TR2 

312  F1=FREQ/FCF1*. 86602540378443900  TR2 

RM1  = RM1 +HA  FI  TR2 

HI =HAF1 *2 , OO  TR2 

GO  TO  ( 3123,  314,3125)  ,NUSE  TR2 

3123  IF (FLAG2)  GO  TO  134  TR2 

FL AG2= . TRUE.  TR2 

GO  TO  1325  TR2 

3125  FLAG3=.TRUE.  TR2 

I NT=0  TR2 

GO  TO  225  TR2 

313  Fl=FREO/FCF2  TR2 

H1=HAF2  TR2 

IF(  NUSE.EQ.3)  GO  TO  3 14  T R2 

IF (FLAG3)  GO  TO  234  TR2 

FLAG3=. TRUE.  TR2 

GO  TO  2325  TR2 

314  I NTF2=0  TR2 

IRCYCL=0  TR2 

315  C0PHC3=PC3( 3) /RM  TR2 

SIPHC3=DSQRT ( 1 . OO-COPHC 3**2)  TRP 

NNUSE(INUSE) =NJSE  TR2 

FCREF  (INUSE)  =FREQ/F  1 TR2 

C CORRECTION  IN  FI  CRITICAL  PRINT  OUT  T R2 

IF ( NUSE . E 0. 2)  FCREF  (INUSE)=FCF1  TR  2 

COT HC3=PC3(l)/0 SORT  (PC3( 1 ) **2  + PC3( 2)  **2)  TR2 

SITHC3=PC  3 (2) /OSQRT  (PC3  ( 1) ** 2+ PC  3 ( 2)  * * 2)  TR2 

PR1 (1)=RTM1*C0THC3*SIPHC3  TR2 
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PR1(2)=RTM1*SITHC3*SIPHC3  TR2 

PR1 (3)=RTM1*C0PHC3  TR2 

Nl(  1)  = COTHC3*SIPHC3»DNOR-COTHC3*COPHC3*ONONL/RTMH-SITHC3*DNOWL/  (RT  TR2 
*M1*SIPHC3)  T R2 

Nl( 2) =S ITHC3*SIPHC3*0NDR-SITHC3*C0PHC3*DN0NL/RTH1-C0TMC3*0N0WL/ (RT  TR2 
»M1*SIPHC3)  TR2 

N1  <3)=COPHC3+ONDR+SIPHC3*ONDNL/RTM1  TR2 

MN1=MAG(N1)  TR2 

T0=-RTM1/HN1  TR2 

1ST  OP  = l TR2 

00  320  J = 1 ,3  TR2 

PO  ( J)  =N1<  Jl  *T0  + PR1(  J)  T R2 

320  A2(J) =P10 (J) -PO (J)  TR2 

CALL  CROSS! A2,C2,XI)  TR2 

CALL  CP0SS(XI,C2jKP)  TR2 

MA2C?  = HAG( XI ) T R2 

HC2  = MAG(C2)  T R2 

00  330  J=l,3  TR2 

M(1,J)=XI(J) /MA2C2  TR2 

H (2,  J)  = C2< J) /MC2  TR2 

330  M(3»J)=KP(J)/(iA2C2*MC2)  TR2 

210PP=M(3,1)*A2(1)+M(3,2)*A2<2)+N(3,  3)*A2(3)  T R2 

RB=  RM1-H 1 TR2 

IF (PB.GE.  OA 8S ( 7 1 OPP)  ) GO  TO  332  TR2 

1ST  0P=2  TR2 

GO  TO  3345  TR2 

332  COBETB=  nABS ( 710PP) /RB  T R2 

SIBFTB=OSQRT ( 1. 00-C08ETB**2>  TR2 

YS1PP=RB*STBETB  TR2 

00  333  J=1  ,3  TR2 

333  S1(J)=M(2, J)*YS1PP  + M(3, J)*Z10PP  + PO(  J)  TR2 

P4=  M AG ( SI ) TR2 

PP2=1.D0- (HR3/R4*C0S02) **2  TR2 

IF(PP2.GE.O.OO)  GO  TO  3331  TR2 

1ST  0P=?  TR2 

GO  TO  3345  T R2 

3331  PP2=R4*0S0RT (PP2) -MR3*SIN02  TR2 

PP2 A=PP?  TR2 

IF(IREFL.GT.l)  GO  TO  3336  TR2 

I F (HITS ) GO  TO  3334  T R2 

PPO=R4*OSQRT  (1.  DO  - (MR3/R4*C0S02)*',2)  TR2 

GO  TO  3336  TR2 

3334  PP0=PP2  TR2 

3336  CBTTM1=RB*C0BETB/RTM1  TR2 

IF( CBTT  Hi . LE • 1 • DO ) GO  TO  334  T R2 

I ST0P=3  TR2 

GO  TO  3345  TR2 

334  SBTTH1=0SQRT  (1.  D0-C9TTM1  + + 2)  T R2 

FLAG=. FALSE.  TR2 

IF( l.OO-Fl+SBTTMl.GT. 0. OO)  GO  TO  335  TR2 

FL AG= .TRUE . TR2 

3345  THETA=1 .DO/RTMJ  TR2 

DIST1  = 1 .010  TR2 

IRCYCL=  0 TR2 

GO  TO  336  TR2 

335  FLOG=H1*F1*OLOG< ( 1.00 +F1*SBTTH1) / ( 1 .00-F1*SBTTM1>) / 2.  OO  TR2 

GPF2=FLOG*2.DO  TR2 
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THETA =CBTTHl*FLOG/RTHi 

PPF2=Hl*SBTTMi+(l  .D0-Fi**(-2> +3BTTHi**2)*FLOG 
IRCYCL= IRC YCL+1 

33  6 YS2PP=YSlPP*OCOS (THETA)  -Z1QPP*DSIN (THETA) 
ZS2PP=ZiOPP*OCOS  (THETA)+YS1PP*DSIN(THETA) 

IF(FLAG  .OR.  I STOP  .GT  .1)  GO  TO  475 

YS3  PP=YS1  PP*  DC  OS  (2.  D0*THETA  ) -ZIOPP’DST^  ( 2 .0  0*THET  A ) 
ZS3PP=Z1OPP*0COS(2 .DO  *THETA)  +YS1PP*DSIN(?.D0*THETA) 
ALPHAV=PI/2.  OQ-THETA-OASIN(SIBETB) 

RV=RB*C0BET8/0SIN ( ALPHA V) 

YVPP=YS2PP*RV/RB 
C32PP= YS3PP- YVPP 
C33PP=ZS3PP-Z10PP 
00  340  J=i»  3 

C3(J)=M(2,J>  *C32PP+M(3,J)*C33PP 
340  A3  (J)=H  (2,J)*YS3PP  + M(3,  J)*ZS3PP*-P0(J) 

CALL  CR0SS(A3,C3 , PHI ) 

CALL  CR0SS(PHI,C3, KPPP) 

M A3C3=MAG ( PH  I) 

MC3=MAG(C3) 

C32PPP=MC3 
DO  350  J=l,3 
P(1,J)=PHI(J)/MA3C3 
P(2,J) =C3 ( J) /MC3 

350  P (3  ,J)=KPPP(J)/  (MA3C3  *MC  3 ) 

C FI  LAYER  OOHNHARD 

YS3PPP=P(2,1)*A3(1) +P  (2, 2) * A3 (2) + P (2, 3)* A3 (3) 

IF(  YS3PPP.GE.0. DO)  GO  TO  3555 
ZS3PPP=P(3,1)*A3(1) *P  (3, 2)*A3(2)+P(3,3)  *A3(3) 

IHT=0 

GO  TO (436 , 3987,353) ,IREFL 
353  IF(MR3.GE.DA9S(ZS3PPP))  GO  TO  360 
IF(B.GE.  OABS ( ZS3PPP) ) GO  TO  356 
3555  IST0P=3 
IRCYCL=0 
GO  TO  475 

356  Y0PPP=-DSQRT(B**2-ZS3PPP**2> 

GO  TO  370 

360  YDPPP=-0SQRT(NR3**2-ZS3PPP**2) 

C FI  LAYER  ITERATION  370  TO  390 
370  00  380  J=1 , 3 

38  0 PO( J)=P(2,J)  *YDPPP+P  (3,J)*ZS3PPP 
NLD=0ASIN(PD(3)/MAG(P0)) 

HL0=HL0N(PD(2)  , PO  (1)) 

CALL  RIIP(RO,NLO,WLD,PLASO> 

RH=HBE+HAE+HAF1*R0 
H=HAF 1 
F=FREQ/FCFi 

IF(RM.  GE.  DABS(ZS3PPP) ) GO  TO  382 
1ST  0P=3 
IRCYCL=  0 
GO  TO  475 

382  Y0PPP2  = -0SQRT(RM**>2-ZS3PPP**2) 

0IST9=0ABS(  Y0PPP2-Y0PPP) 

IF(DISTB,LE,EPS3)  GO  TO  390 

IF (INT.LE. l.OR.OISTB.LT .OISTA)  GO  TO  384 

H=HP 
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F=FP  TR2 

RM=RMP  TR2 

YOPPP2=YDPPP  TR2 

OUT  2=OIST  fl  TR2 

IF(IWRITE.EQ.l)  WRIT E ( 6, 2 4 3) OUT  2 , IHOPS  TR2 

GO  TO  300  T R2 

384  YOPPP=YDPPP2  T R2 

DISTfl=OISTB  TR2 

RMP=RM  TR2 

FP=F  TR2 

HP=H  TR2 

INT=INT+1  TR2 

IF( INT.LE.LIMF1D)  GO  TO  370  TR2 

I STOP=4  TR2 

IRCYCL=  0 TR2 

GO  TO  475  TR2 

C E LAYER  DOWNWARO  TR2 

3 90  YM1 PPP=  Y0PPP2  T R2 

HF 1=H  T R2 

HRM l=OSQRT  (YM1PPP+*  2+7S3PPP++2)  TR2 

DO  395  J= 1 i 3 TR2 

395  RN1  (J)=P(2,  J)  *YH1PPP  + P(3,  J) +ZS3PPP  TR2 

SBETM1=  OABS(  DOT  (RN1  ,C  3)  / < MRH 1* MC 3 ) > T R2 

CBETM1=0SQRT  ( 1 . DO - SBE TM1 * S8ETW1 > TR2 

R4=HAG ( A3)  TR2 

PP2=PP2+R4*0SQRT<1.DO-(MRM1/R4*CBETM1)**2>  -MRM1 +SBET  Ml  TR2 

IF(  F+SBETM1.GT  .1 . DO  ) GO  TO  398  TR2 

1ST  OP=5  TR2 

IRCYCL=  0 TR2 

GO  TO  475  TR2 

398  FLOG=DLOG(  ( F*S3 ETM1 +1 .00)  /<  F * S8ETM1- 1 .D  0) ) TR2 

DGP2P=H+(F/2.D0*FLOG-l.O0/SBETMl)  TR2 

TH12=CBETM1*PGP2P/RM  TR2 

PPEFlO=H*(.5D0»SBETMl+F/4.O0*(l.D0-F  + +C-2)  +C  B ETM 1* * 2 ) *F_  OG-1  . DO/ SB  TR2 
*ETH1)  TR2 

0PP2P=PPEF1D  TR2 

RR2=MRM1»SBETM1  TR2 

RR  2C=HRM1 +C8ET 9 1 TR2 

YM2PPP=YM1PPP  + 0C0S  (TH12)  -ZS3  PPP’DSIN  ( TH12)  TR2 

ZM2PPP=  ZS3PPP+DCOS  C TH  12)  +YM  IPPP'DSIN (TH12)  TR2 

MRM2=OSQRT  (Y  H2  PPP*+  2 ♦ ZM2  PPP* * 2 ) T R2 

YY1PPP=CBETM1+  ( YM2  PPP +C BETH1-ZM2PPP  + S BETM1 ) TR2 

ZY1PPP=CBETH1*(ZM2PPP  + CBETM1  + YM2PPP*SBETM1)  TR2 

C42PPP=  YY1PPP-YM  2PPP  TR2 

C43PPP=ZY1PPP-ZM2PPP  TR2 

MRY1S0=  YY1PPP++2+Z Y1PPP**2  T R2 

MC4SQ=C42PPP**2+C43PPP+*2  TR2 

3984  IF(MR2+*2.GE.MRY1SQ)G0  TO  399  TR2 

IFCB++2.GE.WRY1SQ)  GO  TO  3986  TR2 

I ST0P=6  TR2 

I RC  YCL=  0 TR2 

GO  TO  475  TR; 

3986  TE=-OSORT(  (3 + + 2-MRY  IS  Q)  /MC4SQ)  TRf 

RM=B  T R2 

GO  TO  3995  TR< 

3987  CNST=-00T(A3,C3) /MC3+*2  TR< 

C42PPP=P<2, 1)  + C3(  1)  +P(2,2)  *C3(2)  +P(2 ,3)*C3(3>  TR,' 
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C43PPP=P(3,1)  + C3(  1)  +P (3*  21  *C3(2> +P(3 ,3) *C3 (3)  TR5 

YY1PPP=YS3PPP+CNST* C42PPP  TRJ 

ZY1PPP=ZS3PPP+CNST*C43PPP  TR« 

MRYlSQ=YYlPPP++2+ZY 1PPP++2  TR< 

MC4SQ=MC3**2  T R2 

MRM1=0SQRT  (YS3PPP**  2+  ZS3  PPP  + * 2)  TR< 

CBETH1=DSQRT ( MRY1SQ) /MRH1  TR< 

S8ETH1=DSQRT (1 .00-CBETH1++ 2)  TR< 

GO  TO  3984  TR$ 

399  TE=-DSORT  ( (MR2++2-MRY1SQ ) /NC4SQ ) TR? 

3995  YEPPP=C42PPP*TE+ YY1 PPP  TR^ 

ZEPPP=C43PPP+T£+ZY1PPP  TR^ 

C E LAYER  ITERATION  400  TO  420  TR^ 

INT=0  TR2 

400  OO  410  J = l,3  TR2 

410  PEt  J)=P(2,  J)  *YEPPP+P(3,J)  +ZEPPP  TR2 

NLE=OASIN (PE(3)/HAG (PE) ) TR2 

WLE=HL0N(PE(2) ,PE(1))  TR2 

CALL  RIIP(RO,NLE,WLE,PLASD)  TR2 

NOELAY=. FALSE.  TR2 

IFCFCE.  GT.  0.  DO)  GO  TO  41  05  TR2 

FCE=1 .0-2  TR2 

NOELAY= .TRUE.  TR2 

4105  RM=H0E+HAE+RO  TR2 

H=H  AE  T R2 

F=FREQ/FCE  TR2 

IF  ( RM++  2.  GE.  MRdSQ)  GO  TO  411  TR2 

1ST  OP=6  TR2 

I RC  YC  L = 0 TR2 

GO  TO  4 75  T R2 

411  TE=-0SQRT((RH**2-MRY1SQ) /MC4SQ)  T R2 

YEPPP2= C42PPP*TE+YY1PPP  TR2 

ZEPPP2=C43PPP*TE+ZY 1PPP  TR2 

01 STB=OS3RT  ( (YEPPP-YEPPP  2)  **  2+  (ZEPPP-ZEPPP2)  *+2>  TR2 

IF(OISTB.LE.EPS4)  GO  TO  4.20  T R2 

IFdNT.LE.  1, OR. OIST8.  LT.OISTA)  GO  TO  414  TR2 

H=HP  TR2 

F=  FP  TP2 

RM=RMP  TR2 

YEPPP2=YEPPP  TR2 

ZEPPP2=  ZEPPP  TR2 

OUT  2=0IST  A T R2 

IFdWRITE.EQ.  1)  WRITE (6,1 42) 0UT2,IH0PS  TR2 

GO  TO  4 20  TR2 

414  YEPPP=YEPPP2  TR2 

ZEPPP=ZEPPP2  f R2 

RMP=RH  TR2 

FP=F  TR2 

HP=  H TR2 

D 1ST  A=D 1ST  B TR2 

INT =INT  +1  TR2 

IFdNT.LE. LIMED)  GO  TO  400  TR2 

1ST  OP=7  TR2 

IRCYCL=C  TR2 

GO  TO  475  TR2 

4 20  YM 3PPP=  YEPPP  2 TR2 

ZM3PPP=ZFPPP2  TR2 
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HE  = H TR2 

C ENO  OF  E LAYER  TR2 

MRM3=DSQRT<YM3PPP**2*ZM3PPP**2)  TR2 

C8ETM2=MRM1*C8£TM1/MRH3  TR2 

SBETM2=OSQRT<i.OO-CBETM2*CBET12)  TR2 

IF<IREFL.EQ.  2)  PP2=PP2*MRHl *SBETM1-MRM3*SBETM2  TR2 

IF<F*SBETM2.GT.1.D0>  GO  TO  430  TR2 

ISTOP=8  TR2 

IRCYCL=  0 T R2 

GO  TO  475  TR2 

430  IF( .NOT.NOELAY)  GO  TO  4302  TR2 

FLOG=0. 00  TR2 

OPP1P=O.DO  TR2 

TH34=0.00  T R2 

OGP1P=O.DO  TR2 

YN4PPP=YM3PPP  TR2 

Z H4PPP=ZM3PPP  TR2 

HRM4=MRM3  TR2 

GO  TO  4304  TR2 

4302  FLOG=OLOG<  (F*SBETH2-H.OO) /(F*SBETH2-l.O0> ) T R2 

OGP1P=H*(F/2.DO*FLOG-1. 00/SBETH2)  TR2 

TH34=CBETM2*DGP1P/RM  TR2 

TM4PPP=YH3PPP*0C0S(TH34)  -ZM 3PPP*DL  IN  ( TK  34)  TR2 

ZN4 PPP=ZM3PPP*OCOS(  TH34)  ♦ YH3PPP*DSIM < TH34>  TR2 

MRH4=0SQRT(YM4PPP»*2*ZM4PPP**2>  TR2 

DPP1P=H*<  .50  0*SBETM2«-F/4.00*(1.00-F**  -2)  +CBE1  M2**2)  *FLOG-l.  OO/SBE  TR2 
*TM2)  TR2 

4304  RR1  = HRH4-H  TR2 

IF(  IREFL.LT. 3)  GO  TO  4305  TR2 

IF  ( DABS (RR2C) . GT . RR1 ) GO  TO  431  TR2 

TA  U=RR2 -RR1*DSQRT ( 1 .0  0- ( RR  2C/RR1) **2)  TR2 

PPEF10=PPEF10»0PP1P*TAU  TR2 

GPEF1D=OGP1P*OGP2P*TAU  TR2 

4305  YY2PPP=CBETM2*  (YM4PPP*CBETH2-ZM4PPP*S8ETK2>  TR2 

ZY2PPP=CBETM2*(  ZM4PPP*CBETH2*YM4PPP*SBETM2)  TR2 

C52PPP=YY2PPP-YM4PPP  TR2 

C53PPP=  Z Y2PPP-Z  M4  PPP  T R2 

MC5SQ=C52PPP**2+C53PPP**2  TRJ 

SBM4H2=l.O0-(CBETM2/( 1.O0-HE/NRM4) J **2  TRZ 

IF (SBM4H2.GT .0.0  0)  GO  TO  4307  TRJ 

ISTOP=15  TRJ 

IRCYCL=  0 TRc 

GO  TO  475  TRJ 

4307  TAU1ED=HRH4*SBETM2-(HRH4-HE)  *OSQRT  (SBM4H2)  TRJ 

RY=OSQRT (MRY ISO ) TRI 

GO  TO  437  TRI 

431  1ST  0P=9  T RI 

IRCYCL=  o tr; 

GO  TO  475  TR; 

436  H81  = MAG  (A3)  TR; 

SBET1=-00TCA3,C3) 7(MB1*HC3)  TR, 

C8ET1=0SQRT(1  ,00-SBET  1**  2)  TR 

RY=  MBl^CBETl  TR 

MRY1SQ=RY**2  TR 

CNST=HBl/MC3*SbETl  TR 

C52PPP=P(2,1)  *C3 ( 1) *P(2,2>»C3<2»*P(2  ,3)*C3 (31  TR 

C53PPP=P(3,1)*C3<1)  +P(3,  2)  *C3t  21  *P<  3,  3>  *C3(  31  TR 
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437 

438 

440 

450 

460 
C 

461 

462 

463 
470 

473 

475 

480 

490 

4916 

4917 


Y Y 2PPP- YS3PPP»-CNST*C5  2PPP 

ZY2PPP=7S3PPP+CNST*C53PPP 

MC  5 SQ  = MC3  **Z 

HITS=RO.GT.RY 

I NCHOP=  < RO  +50  .00)  .GE.RY 

IF<RO+HBE.GT„ RY)  GO  TO  438 

I STOP=l 5 

IRCYCL=  0 

GO  TO  475 

TQ=-DSQRT(<(RO«-HBE)  **2-MRYlSQ>  /HC5SQ) 

YBEPPP=C52PPP*TO+YY2PPP 
ZBEPPP=C5  3PPP*T04-ZY2PPP 
IF(  .NOT. HITS)  GO  TO  450 
TQ=-DSQPT( <RO**2-MRY1SO) /HC5S0) 

YOPPP=C52PPP*T9fYY2PPP 
ZQPPP=C53PPP*TQ+Z Y2PPP 
DO  440  J=i  ,3 

PQ(J)=P<2,  J) *YQPPP+ P<3,  J) *ZOPPP 
GO  TO  461 
DO  460  J=l,3 

PQ(  J)  =P(2,  J)  *YY2PPP+P  (3,  J)  *ZY2PPP 
STARTS  TO  ITERATE  TILTED  F2 

0  1ST  2=0 SORT  ( ( PQ(1)  -PREV(l)  ) **2  *■  ( PQ  (2  ) -PREV  ( 2)  >**  2+  (PO  ( 3)  -PREV 
M3)  )**2) 

RT0IF=DABS(RT1-PREV<  13)  ) 

IF ( IRCYCL*  LE. 2)  GO  TO  475 
IF(0IST2-EPS5> 463,463,462 
IF (OIST 1-DIST2) 470 ,475,475 
IF  (RTDIF-FPS5)  500,  500,475 
CALL  PRVGET 
ICNPRB= 1 

IE(IWRI TE.EQ.  1)  WRIT  E(  6,  473)  IHOPS ,E PS5 , DI ST1 , RTD IF 

FORMAT ( 48H  TILTEO  LAYER  CONVERGENCE  PROBLEMS  ON  HOP  NUMBER, 13, 

1 1H,,F6.1,29H  KM.  CRITERIA  NOT  MET,  0IST1=,F8. 2, 8H,  RTDIF= , 

2 F8.2) 

GO  TO  500 

YW2PP=RTM1*YS2PP/RB 
ZW?PP=RTM1*ZS2PP/RB 
00  480  J= 1 , 3 

PH  2 (J)  = M(2,J )*YW2PP+M (3 , J > * ZW2PP+P0 (J) 

NLH2=DASIN(*>W2(3)  /MAG(PW2>) 

WLW  2=  WL  ON  ( PW  2 ( 2 ) ,PW2(1)  ) 

DO  490  J=1 , 3 
W2(J)=PW2(J)-P0 (J) 

BETA22=OACOS(OOT (PW2,  W2)/(MAG(W2) *MAG(PH2>> ) 

BET  TM1=  OASIN ( SBTTM1 ) 

CALL  PRVSTO 
OIST1=OIST2 
BP  = RT  Ml 

CALL  HT3ROF(NOFLAG,FREQ,RTD*NLW2vRTO*HLW2,RTO*BETTM1,RTD*BETA22,B 
*P , RM1 , RTl  , RT  Ml , 0 NOR  , ONONL  , ON  OWL  ,?>N  TFL  G) 

IF< IREFL.NE.  NUSE)  GO  TO  501 
GO  TO(4916, 4917, 4918) ,IREFL 
Fl=  FREO/FCE 
H1=HAE 
GO  TO  4919 

F1=FREQ/FCF1*. 866 0254 037844 3900 
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RM1=RH1*HAF1 
Hl=HflFl*2.00 
GO  TO  4 919 

4918  F1=FREQ/FCF2 
Hl=  HAF2 

4919  INTF2=INTF2+1 
RM=MAG{ PH2) 

00  492  J=l,3 

492  PC3 ( J)=PW2( J) 

IF (INTF2.LE.LIHF2)  GO  TO  315 
IF(ISTOP.GT.l)  GO  TO  7310 
494  IF(FLAG.OR,PNTFLG.GT.5O0)  GO  TO  501 
GO  TO  7310 

C ITERATION  FINISHEO 
500  DO  5003  J = l»  3 
VEC(J,11 =L(2, Jt 

LOC<J,i)=L(2,J)*YBEP*L<3 ,J)*P1PC3> 

VEC  (J,2)=L(2,J)*C12P*L(3,  J) *C13P 
LOC  (J»2  ) = R2(  J) 

VEC  (J»3)=M(2  j J) 

LOCC J,3)=L(2,J)  *Y9P+L<3,J)*Z9P 
VEC  (J,4)=P(2,J) 

LOC  I J»41 =M (2t  J)  *YS1PP*M(3,J)*Z10PP*P0  (J ) 

LOC  (J  ,5  > = A3  (J  ) 

LOC(  J,6i  =P(2  , J)  *YH1PPP+P(3,J)  *ZS3PPP 
•'ECCJ,5)=P(2,J)*C42PPP+P(3,J>*C4  3PPP 
LOC ( Jt  T I =P ( 2 * J)  *YM3PPP*P(3,J) *ZM3PPP 
VEC  (J,E)=P(2,J)*C52PPP*P (3,J)*C53PPP 
50  0 3 LOC  C J»8 )=  P(2 , J)  *YBEPPP*P(  3,  J)  *ZBEPPP 
00  5008  K=1 , 6 
M AG VEC=  HAG ( VEC( 1 1 K) ) 

IF(MAGVEC.EQ.O.DO)  GO  TO  5008 
00  5007  J=l,3 

500  7 VEC<  J,K)  = VFC<  J,K> /MAGVEC 
5008  CONTINUE 

IF (PNTFLG.LT . . 5)  GO  TO  510 

501  IF (FLAG2)  GO  TO  134 
IF!  FLAG  3)  GO  TO  234 
RM=RM1 

INT  = 0 

FLAG=.TRUE. 

GO  T0(502,503,504),  IREFL 

502  C8ET l=COBETB*RB/RN 
FL AG2= • TRUE. 

SBET 1=0  SORT (1.  OO-CBET 1**2 ) 

Y3P=OSQRT(RM**2-Z10PP**2) 

PI  P (3 ) =Z1 0 PP 
F=FREQ/FCE 
H=H  AE 
GO  TO  131 

503  F=  FREQ/FCF1 
FL AG3=. TRUE. 

H=H  AF 1 
MR3=RH 

CBET4=COBETB*RB/RM 

IFCOABS(CBET4)  .LT.l  .OO.ANO.DABS (RM> . GT . DA BS ( Z10PP) ) GOTO  5035 
IST0P=6 
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IRC YCL=  0 TR 

GO  TO  475  TR 

5035  SBET4=OSQRT(1.DO-C8ET4**2)  TR 

Y7P=OSQRT(RM**2-Z10PP**2)  TR 

Z7P=Z10PP  TR 

GO  TO  231  TR 

504  IF(IHRITE.EQ.1)WRITE(6,5Q5)IH0PS  TR 

505  FORMAT ( 23H  PENETRATION  ON  HOP  N0.,I3,1H.)  TR 

GO  TO  7311  TR 

510  IF  CHITS ) GO  TO  530  TR 

GO  TO  (512,514,514) , IREFL  TR. 

512  PP0P=MB1*S8ET1  T R; 

GO  TO  516  TR: 

514  PPOP=(MRM4-HE)*OSQRT(l.DO-(CBETM2/(i.DO-HE/HRM4>)**2)  TR. 

516  DO  520  J-1,3  TR; 

Bl< J)=P (2 , J) *C52PPP  *P (3,J)*C53 PPP  TR! 

520  A1(J)  = PQ(J>  TRJ 

GO  TO  550  TR! 

530  GO  TO  <531, 532, 532), IREFL  TR« 

531  ALPH AQ=OA  SIN (MB1/R0  *C  BET1)  TRZ 

PP0P=MBl*SBETl-R0*DSQRT<l.O0-(RY/R0) **2>  TR1 

GO  TO  533  T Re 

532  ALPHAQ=OASIN(HRM1*CBETM1/RO)  TR2 

PP0P= (MRH4-HEI *OSQRT( 1.0 0- (C8ETM2/ ( 1 . 00-HE/MRH4 ) ) **2 )- RO *DSQ RT ( 1 . 0 TR2 

*0-(MRH4/R0*C8ETM2)**2)  T R2 

533  B22PPP=-C52PPP*OCOS(2.DO*ALPHAQ) +C53 PPP*DSIN( 2. OO’ALPHAO)  TR2 

B23PPP=-C53PPP*OCOS(2.O0*ALPHAQ)-C52PPP*DSIN(2.D0*ALPHAQ)  TR2 

00  540  J= 1,3  TR2 

Al(J)=PO(J)  TR2 

540  B1  (J)=P  (2  ,J)*B22PPP*P  (3,  J)*B23PPP  T R2 

550  IF< IHOPS.GE. IABS  (NODE) ) GO  TO  552  T R2 

551  GO  TO  (5511, 5512, 5513), IREFL  TR2 

5511  PPTOT=PPTOT+PPQ*PPF2+PPOP  TR2 

GPTOT=GPTOT*PPO+GPF2»PPOP  TR2 

GO  TO  5514  TR2 

5512  PPTOT=PPTOT*PPO+DPP1+PPF2+PPZ+DPP1P+PPOP+TAU1EA*TAU1ED  TR2 

GPTOT=GPTOT+PPO  + OGP1«-GPF2*PP2*DGP1P«-PPOPHAU1EA»TAU1EO  T R2 

GO  TO  5 514  TR2 

5513  PPTOT  = PPT OT* PP04PPEF1  A* PPF2 +PP2 +PPEF1 D+PPOP  TR2 

GPTOT=GPTOT+PPO*GPEF1A+GPF2*PP24-GPEF1D4-  PPOP  TR2 

5514  MB1=MAG (Al)  TR2 

IF( .NOT. HITS. ANO, IHOPS.GE. IABS(MODE)  .AND. MB1.GE.RTARG)  GO  TO  598  TR2 
GO  TO  90  TR2 

55  2 IF (SEARCH. AND. RT ARG .GT. MR3-HF1 A . A NO. RTARG ,LE . DMA  XI (RV , MR  3)  ) GOTO  595  TR2 

IF(ASC)  GO  TO  566  TR2 

IF(OABS (RTARG-RO).LE. 1.0-5. ANO. HITS)  GO  TO  580  T R2 

IFC  MOOE.GT . 0)  GO  TO  563  TR2 

C DESCENDING  MODE  TR2 

GO  TO(564, 562,  560)  , IREFL  TR2 

560  IF (RTARG. LE.OMAX1 (RV, MRM1 ) .AND. RTARG. GT .HRMi-HFl)  GO  TO  576  TR2 

IF (DHIN 1 ( MRM4-HE , MRM1 ) . L T . RTARG . AND . RTARG . LE . OMAX 1 ( MRH 4, MRM 1-HF1 ) ) TR2 
* GO  TO  573  TR2 

IF (RTARG.LE.HRN4“HE)  GO  TO  568  TR2 

IF(  INCHOP)  GO  TO  566  T R2 

SEARCH=,TRUE.  TR2 

GO  TO  5513  TR2 

562  IF(RTARG.LE.RV. ANO. RTARG, GT.MRH4)  GO  TO  576  TR2 
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IF(RTARG.L£.HRH4)  GO  TO  568  TR2 

IF( INCH  OP)  GO  TO  566  TR2 

SEARCH= .TRUE.  TR2 

GO  TO  5512  TR2 

563  SEARCH= .TRUE.  TR2 

ASC=. TRUE.  TR2 

GO  TO  551  TR2 

564  IF(RTARG.LE.RV)  GO  TO  568  TR2 

SE  ARCH= • T RUE • T R2 

IF(. NOT. INCHOP)  GO  TO  5511  T R2 

566  1ST  OP=l 5 TR2 

GO  TO  7310  TR2 

568  IFIRTARG.GE. RY. OR. INCHOP)  GO  TO  569  TR2 

SEARCH= .TRUE.  T R2 

GO  TO  551  TR2 

569  TTARG=0.00  TR2 

IFtRTARG.GT. RY) TT ARG=-OSQRT { (R TAR G** 2-M RY 1SQ) /HC 5S Q)  T R2 

YTAPGP=C52PPP*TT  ARG  + YY2PPP  TR2 

ZTARGP=C53PPP*TTARG+ZY2PPP  TR2 

DO  570  J=l,3  TR2 

81 ( J)=C52PPP*P(2, J) +C53PPP*P(3,J)  TR2 

570  PTARG(J)  = P(2,J)*YTARGP+P(3,  J)*ZTARGP  TR2 

PP  0P=  0 • D 0 TR2 

IF(RTARG.GT.RY)PPOP=-RTARG*OSQRT (l.DO-tRY/RTARG) **2)  TR2 

GO  TO  (571 ,5715,5715) ,IREFL  TR2 

571  PP0P=MB1*SBET1*PP0P  TR2 

PPTOT=PPTOT*PPO*PPF2*PPOP  TR2 

GPTOT=GPTOT  + PPO+GPF2  + PPOP  T R2 

GO  TO  600  TR2 

57 15  PP0P= (HRM4-HE) *OSQRT<  1.0 0- < CBET M2/ < 1 . DO- HE /MRM4 ) ) **2 ) ♦ PPO P TR2 

GO  TO(5720, 5720, 5725), IREFL  TR2 

5720  PPTOT=PPTOT*PPO*OPP1+PPF2*PP2*DPP1P+PPOP+T AU1EA+TAU1ED  TR2 

G PT  OT  = GPT OT ♦ PPO ♦ OGP1 *GPF2  *PP2 *OGP IP  *PP  OP  *T  AU 1EA  + T AU 1EO  T R2 

go  to  600  tr; 

5725  PPTOT=PPTOT  + PPO  + PPEF1A*PPF2*PP2  + PPEF1D*PPOP  TR,' 

GPTOT=GPTOT+PPO+GPEF1 A*GPF2*PP2*GPEF1D*PP0P  TRI 

GO  TO  600  TR,' 

573  TTARG=-OSQRT(  ( RTARG** 2-MR Yi SQ> /HC4SQ)  T R,' 

YT  ARGP=  C42PPP  *TT  ARG*YY1 PPP  TR 

ZTARGP=C43PPP*TTARG+ZY1PPP  TR. 

DO  575  J= 1 , 3 TR 

B1 (J)=C42PPP*P(2, J) *C43PPP*P(3, J)  TR 

575  PTARG(J)=P(2, J) *YTARGP*P(3,J) *ZTARGP  TR 

PPTOT=PPTOT*PPEF1A*PPF2+PP2*PPO*PPEF1D  TR 

GPTOT=GPTOT*GPEF1A*GPF2*PP2+PPO*GPEF1D  TR 

GO  TO  600  TR 

576  YT ARGP=-DSQRT ( RT ARG **2- Z S3 PPP* *2 > TR 

ZT  ARGP=  ZS3PPP  TR 

00  579  J=l, 3 TR 

B1(J) =C32PPP*P(2,J)  TR 

579  PTARGCJ)  = P(2 , J) *YTAPGP+P (3,  J)*ZTARGP  TR 

GO  T0<5793, 5793, 5797) , IREFL  TR 

5793  TAU=MRM4*SBETH2-(MRM4-HE)*0SQRT  < 1 . DO - (CBETH2/ (1 . 00-HE/ MRM4 ) ) **2) +P  TR 

* PO*  PP2  TF 

PPTOT=PPTOT»T AU  + OPP1*  PPF2*T AU1EA  TF 

GPT0T=GPT0T*TAU+0GP1+GPF2*TAU1EA  TF 

GO  TO  600  TF 
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5797  TAU=MRM1*SBETH1-(MRH1-HF1)+0SQRT(1. 00-<CBETM1/ (1.D0-HF1/HRM4) ) **2>  TF 

*+PP0+PP2  TR 

PPT  0T=PPT0T+TAU  + PPEF1A+PPF2  +QPP  2P  T F 

GPT  OT=G  PT  OT  + T AU+  GPEF1  A+GPF2  + 0GP2P  TF 

GO  TO  600  TR 

C TARGET  ON  GROUNO  TR 

580  00  590  J=l,3  TR 

B1 <J)=C52PPP*P  (2,  J)  +C53PPP*P<  3,  J)  T R 

590  PT  ARG ( J) = A1 ( J ) TR 

GO  TO(5902, 5904, 5906) ,IREFL  TR 

5902  PPTOT=PPTOT+PPO+PPF2'PPOP  TR 

GPTOT=GPTOT  + PPO  + GPF2-'PPOP  TR 

GO  TO  600  TR 

5904  PPTOT=PPTOT+PPO+OPP1+PPF2+PP2+OPP1P+PPOP+TAU1EA+TAU1ED  TR 

GPTOT=GPTOT+PPO  + OGP 1 + GPF 2+PP2+OG P1P+ PPOP  + T AU1E A + T AU1EO  TR; 

GO  TO  600  tr; 

5906  PPTOT=PPT OT+  PP0+  PPE  Fi A+  PP2+  PPF2  +PPEF1 D+PPQP  TR; 

GPTOT=GPTOT+PPO+GPEF1A+PP2+GPF2+GPEF1D+  PPOP  tr; 

GO  TO  600  tr; 

C ASCENDING  MODE  TRt 

591  YTARGP=OSQRT < RT A RG* *2-Pl P(3 >** 2 > TR2 

00  592  J= 1 , 3 TR< 

592  PTARG(J)=L(2,J)*YTARGP+L<3, J)*P1P(3>  TRZ 

T AU=RTARG*OSQRT ( 1.00- (RY/RT ARG)  **2>  TRt 

PPT0T=PPT0T+TAU  TR2 

GPTOT=GPTOT+TAU  TR2 

INUSE=INUSE-1  TR2 

IF( .NOT  .HITS)  GO  TO  600  TR2 

TAU=-RO*OSORT(1.DO-(RY/RO)**2)  TR2 

PPTOT=PPTOT+TAU  TR2 

GPTOT=GPTOT+TAU  T R2 

GO  TO  600  TR2 

593  TT  ARG=OSQRT( ( RT  ARG**2-MRX1SQ)  /MC1PSQ)  TR2 

YT  ARGP  = C12P*  TTARG  + Y6P  TR2 

ZTARGP=C13P*TTARG+Z6P  TR2 

00  594  J=1 , 3 TR2 

PT  ARG ( J) =L ( 2 , J) * YTARGP+L (3 , J)  *ZTA  RGP  TR2 

594  Bl(J)=L<2,J)*C12P+l<3,J)*C13P  TR2 

T AU  = ( MR  3-HF1 A) *OSQRT ( 1.  DO  - < C0S02/ { 1 . OO-HFl A/MR3) ) * * 2>  TR2 

PPTOT=PPTbT  + OP:>l+TAU  TR2 

GPT  OT=G  PT  OT+  OGP1  + TA  U T R2 

INUSE=INUSE-1  TR2 

IF( .NOT .HITS)  GO  TO  600  TR2 

TAU=-RO*OSQRT(l.DO-( (MR2-HEA) /R0*C0S01>**2)  TR2 

PPTOT=PPTOT+TAU  TR2 

GPTOT=GPTOT  + TAU  T R2 

GO  TO  600  TR2 

595  TTARG=OSQRT< ( RT ARG **2-Yl OP** 2-  Z 10P**2> / CC22P**2 +C23 P**2)  ) TR2 

YT  ARGP  = C22P*TT  ARG  + Y10P  TR2 

ZTARGP=C23P*TTARG+Z10P  TR2 

00  596  J=1  ,3  TR2 

PTARG(J)=L<2, J)  *YTARGP+L(3,J)  *ZTA RGP  TR2 

596  B1CJ)=L(2,J)*C12P+L(3,J)*C13P  TR2 

TAU=(MR2-HEA) +0SQRT (l.OO-tCOSOl/ (1 .00-HEA/MR2))**2) +PP2A  TR2 

PPTOT=PPTOT  +PP  EF 1A  +PPF2 /2. 0 0+T  A J TR2 

GPTOT=GPTOT+GPEFlA  + GPF2/2.  OO  + TAU  TR2 

I F ( . NOT . AS C)  GO  TO  600  TR2 
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INUSE  = INUSE- 1 TR2 

T AU=-R0  +OSQRT (1 .0  0- ( ( MR 2 -HE A)  /RO+COSOl) ++2 ) TR2 

PPTOT=PPTOT+TAU  TR2 

GPTOT=GPTOT+TAU  TR2 

GO  TO  600  T R2 

598  RS  2=M  AG  (A1 ) T R 2 

00  599  1=1,3  TR2 

599  PTARG(I)=A1(I)  TR2 

GO  TO  601  T R2 

600  IF(. NOT. HITS)  RS2  = MAG(PTARG)  T R2 

601  IF (GPONLY)  RETURN  TR2 

NLTARG=RTO+OASIN  <PTARG(3> /RS2)  TR2 

WLT  ARG= RTO+MLON(PTARG(2)  , PTARGC 1) > TR2 

EAST  ( 1)  =-PT  ARG<2)  TR2 

EAST(2)=PTARG(1)  TR2 

EAST( 3) =0. DO  TR2 

CALL  CROSS<PTARG, EAST, NORTH)  TR2 

CNSTl=OOT ( PT  A RG , 01) /RS2  TR2 

ELF=RTD+0ASIN(CNST1/MAG(B1)>  TR2 

ME=MAG( EAST)  TR2 

MN=  HAG ( NORTH)  T R2 

CNST=OOT (EAST ,B1)/HE  TR2 

AZF=RTO*DATAN2(CNST,OOT(NORTH,B1)/HN) + 18 0 . OO-DSIGN ( 18  0 .DO  , CN ST)  TR2 
Bl(  1)=PTARG(  i ) — A 1 1 < 1)  TR2 

B1C2) =PTARG(2)-A11(2)  TR2 

01  (3)=PTARG(3)-A11 (3)  T R2 

S RANGE=  HAG (01)  TR2 

ELSR=RTO+OASIN{OOT (A11,B1)  /(SRANGE  + MA1) ) TR2 

DIST=RO*OACOS(OOT(A11,PTARG)/  (RS2  + MA1) ) TR2 

CNST=OOTCEAST1,B1)  TR2 

BEAR=RT0*0ATAN2(CNST,00T(N0RT-(1,B1)  ) +1B0.D0-DSIGN(  18  0.D0,CNST)  T R2 

RETURN  T R2 

7310  IGOBAK= 1 TR< 

IF(IWRITE.EQ.1)WRITE(6,73G9)  MESSPRC ISTOP)  TRf 

7309  FORMAT  ( 11H  MESSUP—  ,A10)  T R< 

RETURN  TRc 

7311  I GOBAK=  2 TRI 

RETURN  TRI 

END  tr; 
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DOUBLE  PRECISION  FUNCTION  MAG  (A)  MAG 

IMPLICIT  DOUBLE  PRE CISION ( A- Z)  MAG 

DIMENSION  A (3)  MAG 

MAG  = OSORT(A  <1 )+A (1) +A  (2)  *A  < 2)  +A  < 3 > * A ( 3)  > MAG 

RETURN  MAG 

END  MAG 
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DOUBLE  PRECISION  FUNCTION  DOT (A, 8) 

DOT 

IMPLICIT  DOUBLE  PRECISION  <A-H,0-Z) 

DOT 

DIMENSION  A ( 3) ,8(3) 

DOT 

00T=A(1)*B(1)+A(2J*B(2)+A(3)+B(3) 

DOT 

RETURN 

DOT 

END 

DOT 
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SUBROUTINE  CROSS(A,B,C> 

IMPLICIT  (DOUBLE  PRECISION  <A-H»0-Z> 
DIMENSION  A C 3 ) ,B<3>  »C<3) 

C (1)=A(2)*B(3)-A(3)*B  (2) 
C(2)=A(3)*B(1)-A(1)*B(3) 

C(3)  = A(1)*B(2)-A(2)*B  (1) 

RETURN 

END 


VEC 

VEC 

VEC 

VEC 

VEC 

VEC 

VEC 

VEC 
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SUB  ROUT INE  PRVSTO 
DOUBLE  PRECISION  P1,P2 
LOGICAL  LP1,LP2 

COMMON/ CP  RE V/Pl (56) ,P2(56) ,LP1 (6> ,LP2 (6) 
DO  i 1=1,56 

1 P2(I)=P1(I) 

DO  2 1=1,6 

2 LP2 (I)=LP1 (I) 

RETURN 

ENTRY  PRVGET 
DO  3 1=1,56 

3 P1<I)=P2<I) 

DO  4 1=1,6 

4 LP1 (I)=LP2 (I) 

RETURN 

END 


ROUTINE  HTGROE  74/74  OPT=l  FTN  4.5+414  04/27/ 


SUBROUTINE  HTGROF (NOFL  AG, FREQ, NL AT , WLONG, ANGLE1 , ANGLE2,  RANG12,  RMU  HTG 
+,RTU,RTMU,ONOR,ONOT,ONCP,PNTFLG>  HTG 

IMPLICIT  DOUBLE  PRECISION  <A-H,0-Z>  HTG 

OOUBLE  PRECISION  NLAT , NLA  TR  HTG 

LOGICAL  TILT, FAST  HTG 

COMMON/ TIL TC /TILT ,FAST  HTG 

NLATR=NLAT/57 .2957795130823200  HTG 

ML ONGP=W LONG/5 7. 2957795 13 08232D0  HTG 

CALL  RIIP( 63 70. DO, NL ATR, WLONG R,EN)  HTG 

CALL  RTFIND(N0FLAG,FREQ,ANGLE1,  ANGLE2,RANG12,RTU,RMU, RTMU,PNTFLG)  HTG 
IF (TILT)  GO  TO  1 HTG 

ONDR=1.000  HTG 

DNOP=  0. OD  0 HTG 

DNOT=  0.  OD 0 HTG 

RETURN  HTG 

1 CALL  OENSEIRTU, NLAT, WLONG, EN,DNDR,DNDT ,DNDP>  HTG 

RETURN  HTG 

END  HTG 


FTN  4.5+414  04/27/ 


PRV 

PRV 

PRV 

PRV 

PRV 

PRV 

PRV 

PRV 

PRV 

PRV 

PRV 

PRV 

PRV 

PRV 

PRV 

PRV 
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IROUTINE  RTFINO 


74/74 


OPT  = i 


FT N 4. 5+414 


04/27/ 


SUBROUTINE  RTFI NO ( NOFLA G,FREQ , AN GL El , AN GLE2 , RANG1 2 , RTU , RMU, RTMU ,P  RTF 


*NTFLG)  RTF 

C NOFLAG  FLAG  TO  INOICATE  WHICH  PASS  RTF 

C =0  FIRST  PASS  RTF 

C =1  SUBSEQUENT  PASSES  RTF 

C NOTE  - THIS  ROUTINE  REQUIRES  INITIALIZATION  . RTF 

C ON  FIRST  PASS,  ANGLE1  MUST  = THE  INITIAL  TAKE  OF-  ANGLE  ON  RTF 

C GROUNO,  ANGLE2  HUST  =0,  RANG12  MUST  = RADIUS  OF  EARTH  RTF 

C FREQ-OPERATING  FREQUENCY  IN  MHZ.  RTF 

C ANGLE1  - ANGLE  (OEG)  BETWEEN  TILTED  LAYER  AND  RAY  INCIDENT  ON  THE  LAY  RTF 

C ANGLE2  - TILT  ANGLE  (DEG)  - THAT  IS,  ANGLE  BETWEEN  TILTED  LAYER  AND  A RTF 

C EXACTLY  HORIZONTAL  LAYER  AT  THE  REFLECTION  POINT  RTF 

C RANG12  - DISTANCE  FROM  CENTER  OF  EARTH  (KM)  OF  REFLECTION  POINT  USED  RTF 

C COMPUTE  ANGLE1  + 2 RTF 

C RTU-COMPUTED  DISTANCE  FROM  CENTER  OF  EARTH  IN  KM.  OF  POINT  OF  INTERES  RTF 

C RMU-DI  ST  AN  CE  FROM  CENTER  OF  EARTH  TO  LAYER  IN  KM.  RTF 

C PNTFLG-PENETRATION  FLAG  RTF 

C =0.  REFLECTION  RTF 

C =1.  PENETRATION  RTF 

IMPLICIT  OOUBLE  PRECISION  (A-H,0-Z)  RTF 

DOUBLE  PRECISION  M3000  RTF 

DIMENSION  RM(3) ,H(3),FC(3) ,RTM(3)  RTF 

COhi’ON/RIIP  AR/M30  0 0 , FCF2  ,FCF1,FCE*HBE,HAE,HME,HAF1,HMF1,HAF2,HMF2*  RTF 
1DUM  (13)  ,10(3)  RTF 

COHMOKY  RTCOM/COSBO,  NUSE  RTF 

DATA  R,\ 0/0. 1745329251 99 4330D-01/  RTF 

D ACOS (DUMMY) =0 AT AN2 (OSQRT ( 1 . 0 OO -DUMMY* DUMMY)  , DUMMY)  RTF 

PNTFLG::0. 0D0  RTF 

IF  (NOFI  AG)  4,  4 , 44  RTF 

4 IF( ANG.E2. EQ.  0. 0 00) GO  TO  44  RTF 

WRITE ( i , 10 0) ANGLE2  RTF 

100  FORMAT  !22H  THE  VALUE  OF  ANGLE2  =,  015.  8,13H  IS  IN  ERROR.)  RTF 

PNTFLG- 1. ODO  RTF 

44  FC(l)=FOE  RTF 

FC(2)=FCC1  RTF 

FC  ( 3)  =FCF  2 RTF 

H(1 )=HAE  RTF 

H(2)  = HAF1  RTF 

H(3)=HAF2  RTF 

RM(3)  =6  37  0 .0  0 +HMF  2 RTF 

RM(2)=RM(3)-H(3)  RTF 

RM ( 1) =RM( 2) -H ( 2)  RTF 

IF ( NOFLAG. GT , 0)  GO  TO  45  RTF 

COSBO= ( OCOS(ANGL El* RAO) *RANG 12/ 6370.00)  RTF 

45  00  10  N=1 , 3 RTF 

C COMPUTE  RTM  FOR  EACH  LAYER  RTF 

IF(FC(N)  .EQ.  0.  OD  O GO  TO  6 RTF 

FAC=(RH  (N  )**2)  /16.D0-  (H ( N)  **  2/ 2.  DO ) *(  ( FREQ/FC  ( N)  ) **2-1 .00)  RTF 

IF(FAC) 6,5,5  RTF 

5 RTM(N)=0.7500*RM(N) +OSQRT (FAC ) RTF 

GO  TO  10  RTF 

6 RTM (N)=RM ( 5)  RTF 

10  CONTINUE  RTF 

C FIND  MAX  RTM  RTF 

RMM  = RTH(3)  RTF 

IF(RTM(3)  . EQ. RM ( 3) ) RMM=RM(3)/2.CD0  RTF 

3000  IF (NOFL AG.GT .0)  GO  TO  14  RTF 
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ROUTINE  RTFINO 


7 4/74 


DPT  = 1 


FT  N 4. 5 + 4L4 


04/27/ 


C IF  ALL  RTM  ARE  EQUAL,  USE  F2  LAYER  RTF 

IF(RTM<2> .NE.RTM(3) .OR.RTM  ( 1) .NE.RTM (3) )GO  TO  12  RTF 

NUSE=  3 RTF 

RTHUSE=RTM( 3)  RTF 

GO  TO  14  RTF 

C FIND  MINIMUM  RTM  RTF 

12  PTMUSE=RTM (1 ) RTF 

NUSE=1  RTF 

00  13  1=2,3  RTF 

IF(RTHUSE.LT.RTMd)  )G0  TO  13  RTF 

RTHUSE= RTM(I)  RTF 

NUSE= I RTF 

13  CONTINUE  RTF 

IF(NDFLAG.GT.O)  GO  TO  14  RTF 

123  ANGLE1= (OACOS(COSBO  *6370  .D0/RANG12) ) /RAD  RTF 

C USE  PARAMETERS  ASSOC.  WITH  MIN  RTM  RTF 

14  RTMUSE=RT M( NUSE)  RTF 

RMUSE=RM(NUSE)  RTF 

HUS  E=H (NUSE)  RTF 

FCUSE=FC(NUSE)  RTF 

144  C8TM=0C0S ( ANG  LEI *RA  0)  *RA  NG12/RTHUSE  RTF 

IF(CBTM-1.  000)  1 5,15  ,16  RTF 

15  0TM=OACOS(C9TM)  RTF 

17  X = 1.0D0-< (FREQ/FCUSE)**2)* < (OSIN(BTM) >**2)  RTF 

IF(X.GE. 0.000) GO  TO  20  RTF 

16  IF(RTMUSE.NE. RM(3) ) GO  TO  18  RTF 

C P ENET  PAT  I ON  RTF 

185  RTU=RMM  RTF 

PNTFLG= 1.000  RTF 

RMU=P.M(  3)  RTF 

RTMU=RTMUSE  RTF 

2000  RETURN  RTF 

18  IF<NUSE-3)19,185,185  RTF 

C ELIMINATE  RTM  USEO  AND  FIND  NEXT  SMALLEST  RTM  RTF 

19  NUS  E=NU  SE* 1 RTF 

GO  TO  123  RTF 

C REFLECTION  RTF 

20  RT=RTHUSE-HUSE*OSCiRT(X)  RTF 

RMU=RMUSE  RTF 

RTMU=  RT  MUSE  RTF 

A = RMU-PT  RTF 

B=  A/DCC  S( ANSL  E2*RAO ) RTF 

IF (B.GE ,HUSE)B=HUSE  RTF 

RTU=RMU-B  RTF 

END  RTF 
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ROUTINE  DENSE  74/74  OPT  = i FTN  4.5+414  04/27/ 

SUBROUTINE  OENS E < R ,T HET A, PHI,E N, ONDR,  DNDT , DNDP)  DEN 

C R-OISTANCE  FROM  CENTER  OF  EARTH  IN  KM.  OEN 

C THETA-NORTH  LATITUOE  IN  DEGREES  DEN 

C PHI-WEST  LONGITUDE  IN  DEGREES  DEN 

C EN-VALUE  OF  ELECTRON  DENSITY  AT  GIVEN  POINT  DEN 

C DNOR-PARTIAL  DERIVATIVE  OF  ELEC.  DEN.  WITH  RESPECT  TO  R OEN 

C ONOT-PARTIAL  DERIVATIVE  OF  ELEC.  DEN.  WITH  RESPECT  TO  THETA  DEN 

C ONDP-PARTIAL  DERIVATIVE  OF  ELEC.  DEN.  WITH  RESPECT  TO  PHI  DEN 

IMPLICIT  OOUBLE  PRECISION  (A-H,0-Z>  DEN 

DIMENSION  R X ( 9 ) ,RY(9)  » ALPHA (9) »FN ( 9)  OEN 

COMMON/RPERT/RRIIP, ITIP  DEN 

DATA  RAD/ . 0174532925199433D0/  OEN 

10  RY(1)=R  DEN 

RX(i)=THETA*RAD  DEN 

ALPHA (1)= PH I+RAO  DEN 

NSIG=1  DEN 

20  CALL  DERV(NSIG,FN,RY,RX, ALPHA, DNDR, DNDT, DNOP)  OEN 

IF (NSIG.EQ .3 ) GO  TO  30  DEN 

DO  25  1=2,5  DEN 

CALL  PTIP (RY (I) ,RX  II) , ALPHA (I) , FN(I)  ) OEN 

GO  TO  (22,23)  ,ITIP  OEN 

22  CALL  NFR0MR(RY(I+4),FN(I+4>>  DEN 

GO  TO  2fc  OEN 

23  CALL  RIIP(RY ( 1 + 4) , RX (1  + 4) , ALPHA (I +4)  , FN(I+4) ) DEN 

25  CONTINUE  OEN 

CALL  RIIP(RY(1) ,RX(1) , ALPHA (1)  ,FN(1))  DEN 

NSIG=2  DEN 

GO  TO  20  OEN 

30  EN=FN(1>  DEN 

RETURN  DEN 

END  DEN 
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DUTINE  DERV 


74/74 


0PT  = 1 


FTN  4.5  + 414 


04/27/ 


SUBROUTINE  OERV( NS IG , FN, RY  , RX , ALPHA , DNOR , DNDT, DNDP) 

C NSIG-SWI TCH  TO  INDICATE  WHETHER  POINTS  OR  DERIVATIVES  ARE  REQUIRED 

C FN-ARRAY  OF  9 VALUES  OF  ELEC.  DEN. 

C RY-ARRAY  OF  9 VALUES  OF  R IN  KM. 

C RX-ARRAY  OF  9 VALUES  OF  THETA  IN  RAD. 

C ALPHA-ARRAY  OF  9 VALUES  OF  PHI  IN  RAD. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  FN(9) ,RY<9) ,RX(9) .ALPHA (9  * ,0F(3) 

DATA  S/1.  0D0/ 

IF  CNSIG-2)  100,  200,200 

C FIND  8 POINTS  ON  CUBIC  SURROUNDING  RY  (1) ,RX(1> , ALPHA (1) 

100  X1=S/(0.17320508D+01) 

X2=X1/RY ( 1) 

C=RY(l)+OCOS(RX(l>) 

IF(C.EQ.0.0O0)C=100.0D0 
X3=X1 /C 
DO  10  1 = 2,5 
RY  (I)=RY(1)+X1 
RY ( 1+4) =RY ( 1) -XI 
10  CONTINUE 

DO  20  1=2,8, 2 
RX(  I)  =PX(  1)  +X2 
RX(I+1)=RX(1)-X2 
20  CONTINUE 

DO  30  1=2,6, 4 
ALPHA (I )=ALPHA ( 1) +X3 
ALPHA ( 1+11  =A LPHA (I) 

ALPHA (I+2)=ALPHA(1)-X3 
ALPHA(I+3)=ALPHA (1+2) 

30  CONTINUE 

RETURN 

C FIND  PARTIALS 
200  AM=0. 125000 

DF(1)=AN*(FN (2)  -FN (61 +FN  ( 3)-FN (7 ) +FN (4)-FN ( 8 ) + FN (5 ) - FN (9> ) 

OF  (2)  = A M+  ( FN  (2)  -FN(3)  +FN(4)  -FN  ( 5)  +FN  ( 6)  -FN  (7)  +FN  ( 8)  - FN  (9)  > 
OF(3)=AM*(FN(2)-FN(4) +FN ( 3) -FN ( 5 ) +FN (6 ) -FN ( 8) +FN(7) -FN(9>> 

240  DNOR=OF(l)/Xl 

DNDT=OF(2)/X2 
ONOP=OF (3)/X3 
NSIG=3 
RETURN 
END 
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Appendix  C 


Block  Diagram,  Flow  Chart,  and  Program  Listing 
for  WIMP  Driving  Program  OBLFACT 


Explanation: 

Cl.  OIU.FACT  BLOCK  DIAGRAM 

This  presentation  illustrates  the  communication  and  iteration  network  employed 

to  generate  the  IU(3000)F,-,  and  f Fg  gradient  correction  factors  required  to  match 

the  predicted  leading  edge  to  the  given  oblique  ionogram.  The  function  of  the  various 

blocks  is  generally  described  as  follows; 

BLOCK  I:  Initialization  and  Input  Procedures. 

BLOCK  11:  Refinement  of  IVRSOOOlFg  Gradient  Correction  Factor. 

BLOCKS  III  and  IV : Refinement  of  f F_  Gradient  Correction  Factor. 

o 2 

BLOCK  V:  Termination  and  Output  Procedures. 

C.2.  DKTAILKI)  FLOW  CHARTS 

Detailed  flow  charts  are  presented  for  each  block.  Symbols  and  conventions 
are  the  same  as  described  in  Appendix  A. 


c;l  program  listing 

Listings  arc  included  for  OBLFACT,  PEAK,  and  SKIP. 
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Block  I,  Detailed  Flowchart 
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Block  III,  Detailed  Flowchart 
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Block  IV,  Detailed  Flowchart 


\L87, 
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(OUTINE  OBL 


74/74 


OP  T =1 


FT  N 4.5  + 414  04/27/ 


SUBROUTINE  OBL 

C PROGRAM  OBLFACT (INPUT, OUTPUT, TAPE5=INPUT,TAPE6=OUTPUT) 

IMPLICIT  OOU8LE  PRECI SI  ON { A -H  ,M ,0 -Z) 

INTEGER  MOOE 

OOUBLE  PRECISION  NLi,NLREF 
LOGICAL  SMOOTH, GPONLY 
COHHON/ GPFLAG/GPONLY, IWRITE 

COMMON/ FIDOL E /FP, FS ,R1,NL1 ,HL1 , A Z,RS, EL  STEP,  MOOE, NSKPCL 
COMMON/ PERTC/NLREF,  WL RE F , PERTP ( 14) 

COMMON /RIIPAR/OUM 1(11) , SN , ZT , YRMO , OOY , H AESET , MF AC, F2FAC, FI  FAC, 
1 EFAC,H2FAC,0UM2(T) ,IO,IA, SMOOTH 
NAHFL IS T/ INPUT/ SN,ZT, YRMO, OOY, MFAC,F2FAC, 

1NL 1, WL 1, Rl, NLREF, WLPEF.RS ,AZ,ELSTEP,  MODE, 
1FCF2R,HMINR,0CHK,FP,0ELP ,FS, DELS, HRO , FRO 
3 , IRSTRT , ZTREF 
IRSTRT  = 1 
GPONLY= . TRUE • 

IWRITE=0 
10=1 
I A=  0 

Rl=  63  70  .00 
RS=  6370 • DO 
MOOE=-1 
NL1  = 

HL1  = 

OCH  K= 10.00 
SN3  0=  0. OO 
NLREF=NL1 
HLREF=HL1 
AZ= 

ELSTEP=  1.00 
1 FR0=0 . 0 0 

MRO  =0 .D  0 


4 


CCCC 

900 


901 


902 


FP=  0.00 
OELP=  0. DO 
REA  0( 5, INPUT) 

IF(IRSTRT.EQ.O)  STOP 
IF(SN30.EQ.YRMO)  GO  TO  4 
SN30= YRMO 

CALL  RIIP(0. 00,0. DO, 0.00, 0.00) 

IPKCNT=  0 

FCF2R=-  DABS (FCF2R) 

HMINR=-  DABS  (HMINR) 

CALL  FACTO (NLREF, WLREF,FCF2R, HMINR, Z TREF) 
HRITE(6, INPUT) 

WRITE (6,90  0)  SN,ZT,VRMO,OOY,MFAC,F2FAC 
FORMAT( IX,  (SN  =( , 11 X, D24 . 18/1X , [ ZT 

* IX, (YRMO  =[,11X, D24.18/1X,  COOY 

* IX , [ MF  AC  =( ,11X,D24.18/1X,IF2FAC 

WRITE(6,901)NL1,WL1,R1,NLREF,WLREF,RS 
FORMAT ( IX , t NL1  =C , 11 X , D24. 18/1X , [ WL1 

* IX, [R1  =C ,11X, 024. 1B/1X, [NLREF 

* IX, [WLREF  =[ ,11X,024.18/1X,[RS 

WRITE (6,90  2)  AZ,  EL  STEP, MOOE, FCF2R, HMINR, DCHK 


FORMAT ( IX , [A  Z 
1 IX, [MOOE 

‘ IX, [HMINR 


=[ ,11X,024.18/1X,[ ELSTEP 
=[  ,1 6X , 13  /IX , [ FCF2R 
=( ,11X,D24.18/1X, (OCHK 


=[,11X,D?4.18/ 
=[, 11X, 024.18/ 
= [,  11X,  D2  4.  18) 

=t , 1 IX, 024. 18/ 
=[, 11X.D24.18/ 
= [,  11X, 024,  18) 

=[  , 11X,  D24.  18/ 
=[ ,11X ,024. 18/ 
=[ ,11X,024, 18) 
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OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 
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OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 
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i ROUT  r N E OBL  74/74  OPT  = i FT  N 4.5*414  04/27 

WRITE(6,903)FP,DELP,  FS.OELS,  MRO , FRO  OBL 

90  3 FORHATdX,  (FP  = [ , 11X, 0 24. 18/ IX , [ DELP  =t ,11X ,024. 18/  OBL 

* IX, (FS  =[  ,11X, 024. 18/1  X, t DELS  =t , 11 X, 024. 18 / OBL 

* IX , [ MRO  =[  ,11X,024.18/1X, [FRO  =[ , 11X, 024. 19)  OBL 

WRITE(6 ,904) IRSTRT, ZTREF  OBL 

904  FORMATdX,  [IRSTRT  =t,16X,I3  /IX,  [ZTREF  =( , 11X,024.  18)  OBL 

IHO  1=  0 OBL 

IFO  1=  0 OBL 

NS  KPOL=  0 OBL 

IGOBAK=  0 OBL 

GPP  0=  0. OO  OBL 

IF(DA9S  (DELP) .LT.DCHK. AND. DABS(OELS)  .LT.DCHK)  GOTO  330  OBL 

IF( FP.GE.5.740Q)  CALL  PEA K (FRO , MRO , GPP  0 ,IGOBAK)  OBL 

IF (IGOBAK. GT. 0)  GO  TO  310  OBL 

CALL  SKIP (FR  0, MR  0 , GPS  0,  0,  IGOBAK)  OBL 

IF ( IGOBAK , GT . 0 ) GO  TO  315  03L 

GPPOES=GPPO*DELP  OBL 

GPSDES=GPSO*OELS  OBL 

WRI  TE  (6 , 1 1 1)  GPPOES , G PSD  ES  OBL 

111  FORMA  T< ( OESIREO  G PP = t , 4P0 1 5. 2 » t » GPS=[,015.2)  OBL 

DELP1=DELP  OBL 

DEL  Sl  = OEL  S OBL 

DELPO=OELP  OBL 

OELS  0=0  ELS  OBL 

MR1=MR0  OBL 

IF(OABS(OELP)  .LT.DCHK)  GO  TO  197  OBL 

GO  TO  107  OBL 

C MR  FIDDLING  OBL 

100  IF(FP.LT. 5.7400)  GO  TO  300  OBL 

CALL  PEAK(FR0, MRO, GPPO, IGOBAK)  OBL 

IF( IGOPAK.E9.0)  GO  TO  105  OBL 

MR l=MRO  OBL 

MR0=MR1 00100  OBL 

IPKCNT=IPKCNTM  OBL 

I F( IPKC  NT-10 0)  100,320,320  OBL 

105  DELPO=GPPOES-GPPO  OBL 

NSTMT  = 105  OBL 

WRITE(6,llll)  NSTMT , FRO, MRO, OELPO  OBL 

1111  FORMA T(  I10,lP3O15 .7)  OBL 

IF(DABS (OELPO)  . LT.  OCHK)  GO  TO  300  OBL 

107  MRl=MR0-DSIGN(.001O0, OELPO)  OBL 

110  CALL  PEAK(FR0, MRI, GPP1, IGOBAK)  OBL 

IFUGOBAK.EQ.  0)  GO  TO  115  OBL 

IF(OELPO.GT.O.OO)  GO  TO  128  OBL 

MR0=HR1  OBL 

MR1=MRO-OSIGN(  .00100, OELPO)  OBL 

OELP0=-10000 .00  OBL 

IPKCNT= IPKCNT  + 1 OBL 

IF(IPKCNT-IOO)  110,320,320  OBL 

115  0ELP1= GPPOES -GPP  1 OBL 

NSTMT=1 15  OBL 

WRITE  (6  ,1111 ) NSTMT, FRO, MRI, OELP1  OBL 

IF(DABSIOELPI)  .LT.OCHK)  GO  TO  190  OBL 

IF ( OELP l*OELP  0 ) 130,190,120  OBL 

120  MRR=DEL  PI* (MRI- MRO)  / ( OELPO- OELP1)  OBL 

IFtOABS(MRR) .GT. 0.100)  MRR=OSI GN ( 0 . ID 0, MRR)  OBL 

MRR=MR1 +MRR  OBL 
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FTN  4.5+414 


04/27/ 


IF(-OELPO.EQ. 1000  0. 00  0.  OR. (HR1-MRR) *OELP0 .L E . 0 .0 0) 

OBL 

1 HRR=HR1-DSIGN(.  00100  ,OELPi) 

OBL 

* 

MR  0=HR1 

OBL 

HR1=MRR 

OBL 

u 

i 

DELPO=OEL  PI 

OBL 

GO  TO  110 

OBL 

128 

O EL Pi =-10000. 00 

OBL 

1 

NSTMT=128 

OBL 

WRITE  (6  ,1111)  NSTMT,FR0,MR1,DELP1 

OBL 

130 

IF(OABS(MR1-MRO>  . LT.  . 000100  ) GO  TO  180 

OBL 

1 

MRR= (MR1+MR0 ) * • 50  0 

OBL 

CALL  PEAK(FR0,MRR,GPPR, IGOBAK) 

OBL 

IF (IGOBAK .GT.O)  GO  TO  160 

OBL 

DELPR=G  PPOES-G  PPR 

OBL 

NSTMT=1 30 

OBL 

WRITE (6, 1111)  NSTMT ,FR0 , MRR,OELPR 

OBL 

IF( DA BS (DELPR) .LT.OCHK)  GO  TO  195 

OBL 

IF(DELPR*OELPO)  150,195,140 

OBL 

140 

IFIDABS (DELPR) .GT.DABS(OELPO) ) GO  TO  180 

OBL 

OELPO=OELPR 

OBL 

1 

MR  0=MRR 

OBL 

GOTO  130 

OBL 

, 

150 

IF< DABS (DELPR)  .GT.OABS(OELPll)  GO  TO  180 

OBL 

OELPl=OELPR 

OBL 

MR1=MRR 

OBL 

GOTO  130 

OBL 

160 

DELP1=-10000.00 

OBL 

HR1=MRR 

OBL 

GO  TO  130 

OBL 

180 

IF(IM0i*IF01.EQ.l)  GOTO  325 

OBL 

IHO 1=1 

OBL 

IF( DABS(OELPl)  .GE.DABS(OELPO) ) GO  TO  200 

OBL 

DELP0=DELP1 

OBL 

MR0=HR1 

OBL 

GOTO  200 

OBL 

190 

MR0=MR1 

OBL 

i 

OELPO  =0ELP1 

OBL 

GOTO  197 

OBL 

! 

195 

MR0=MRR 

OBL 

DELPO=OELPR 

OBL 

[ 

197 

IMO 1=0 

OBL 

1 

C FR 

FIDDLE 

OBL 

200 

CALL  SKIP(FRO,MRO,GPSO,1 , IGOBAK) 

OBL 

l 

IF  (IGOBAK -1)  205,201,  320 

OBL 

201 

FR1=FR0 

OBL 

FR0=FRl+.001O0 

OBL 

i 

GO  TO  200 

OBL 

4 

205 

DELSO=GPSOES-GPS0 

OBL 

NSTMT=205 

OBL 

WRITE (6,1 111)  NSTMT, FRO, MRO,OEL SO 

OBL 

IF(OABS (OELSO). LT.OCHK)  GO  TO  300 

OBL 

• 

FR1=FRO-OSIGN ( . 0 0 10 0 , DELS  0) 

OBL 

210 

CALL  SKIP(FR1,MR0,GPS1,1, IGOBAK) 

OBL 

IF( IGOBAK-1)  215,211,320 

OBL 

211 

IF(DELSO.GT.O.OO)  GO  TO  226 

OBL 

FR0=FR1 

OBL 

\ 

FR1=FRO-OSIGN( .00100, DEL  SO) 

OBL 

L 
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ROUTINE  OBL 


74/74  OPT=l 


FT N 4.5  + 414 


04/27 


OEL  S0  = -10000  .00 
GO  TO  2 10 

215  OELSl=GPSOES-GPSl 
NSTMT=2 15 

WRITE (6, 1111)  NSTHT ,FR1 , MRO  ,0ELS1 
IF(OABS(OELSl) .LT.OCHK)  GO  TO  290 
IF(D£LS1+OELSO>  230,290,220 
220  FRR=D£LS1* (FR1-FR0) / (OELSO-OELS1) 

IF( DABS(FRR)  -GT.O.IDO)  FRR=DSIGN (0. 100  ,FRR) 

FRR=FR1 +FRR 

IFC-OELSO.EO. 10000. 000. OR. (FR1-FRR) *DELS0 .LE.O.OO) 
1 FRR=FR1-DSIGN(,  00100, OELS1) 

FR0=FR1 
FR1=FRR 
OELSO=DELS1 
GOTO  210 

226  QELS1=-10000.00 

NSTMT  = 2 28 

WRITE (6,1111)  NSTHT, F Rl, MR 0,OELSi 
230  IF(OABS(FR1-FRO)  . LT . . 000100)  GO  TO  280 
FRR= (FRl+FRO)’ .500 
CALL  SKTP(FRR,MRO,GPSR,  1, IGOBAK) 

IF(IGOBAK-l)  235,260,320 
235  OELSR=GPSOES- GPS R 
NSTMT=235 

WRITE (6, 1111)  NSTHT, FRR,  HRO, OELSR 


IF(OABS (OELSR) .LT.OCHKI  GO  TO  295 
IF (OEL SR* OEL SO)  250,295,240 

240  IF(OABS  (OELSR).  GT . DAB S(  OEL SO) ) GO  TO  28  0 
OELSO=PELSR 
FRO  =FRR 
GOTO  230 

250  IF (DABS  (OELSR) .GT .DAB  S(OELSl) > GO  TO  280 
OELSl=OELSR 
FR1=FRR 
GOTO  230 

260  OELS1=-10000.00 
FR 1=FRR 
GO  TO  230 

280  IF (IH01+IF01 . EQ . 1)  GO  TO  325 
IFO 1=1 

IF (PASS  (OELSl) .GE.DABS(OELSO) ) GO  TO  100 
DELS  0=0  ELS  1 
FR  Q=FR1 
GOTO  100 

290  FR0=FR1 

OELSO=OELS1 
GOTO  297 

295  FR0=FRR 

OELSO=OELSR 

297  IFO 1=0 
GOTO  100 

300  WRITE  (6,305)  HRO  , FRO  , OE  LPO  , OE  LSO 

305  FORMAT ( t 0 CONVERGENT  VALUESV  HR=[ , 1P015. 7 , 1 , FR=t,015.7/ 
1 1 0 X, t MISS  GIST • PEAK= 1,010.3,1,  S KIP= I , DIO .3) 


OBL 

OBL  ^ 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL 

OBL  t 

OBL  4 

OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
OBL 
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ROUTINE  081 


7 4/  74  0PT  = 1 


FTN  4.5+414 


04/27/ 


GOTO  1 OBL 

310  WRITE (6*312)  OBL 

312  FORMA T ( [ 0 PEAK  RAYTRACE  FAILS  ON  INPUT  CASE,  TRY  DIFFERENT  INPUT!)  OBL 

GOTO  1 OBL 

315  WRITE(6,317)  OBL 

317  F ORMA  T ( 1 0 SKIP  RAYTRACES  FAIL  ON  INPUT  CASE,  TRY  DIFFERENT  INPUT!)  OBL 

GOTO  1 OBL 

320  WRITE  (6  ,322)  OBL 

322  FORMAT! !0  100  ITERATIONS  DONE  BUT  DOESNt , 1HI ,t  T CONVERGE,!,  OBL 

It  TRY  DIFFERENT  INPUT!)  OBL 

GOTO  300  OBL 

325  WRITE(6, 327)  OBL 

327  FORMAT! tO  CONVERGENCE  PROBLEMS,  T RY  DIFFERENT  INPUT!)  OBL 

GOTO  300  OBL 

330  WRITE !6 ,332)  OBL 

332  FORMAT!!  YOUR  ORIGINAL  INPUT  CONVERGES!)  OBL 

GOTO  1 OBL 

END  OBL 


OUTINE  SKIP 


74/74  0PT= 1 FTN  4.5+414  04/27/ 

SUBROUTINE  S KIP ( FRR,MRR , GPS , I OPT , IGO BAK  ) SKI 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z)  SKI 

REAL  OUT ( 3)  SKI 

DOUBLE  PRECISION  MRR, HR, NLi , NLT, NLREF , GP( 42) , RS2 ( 42)  SKI 

COMMON/ PERTC/NLRE F, WLREF, FR , FD,HR , DUM (11)  SKI 

COMMON/FIODLE/FP ,FS ,R1 , NLI , HLi , AZ  , RS, ELSTEP, MOOR, NSKPCL  S KI 

COMMON/ R£MWEN/INUSE,NNJS  E< 20) , FCREF ( 20)  SKI 

NSKPCL=  NSKPCL+1  SKI 

IF (NSKPCL. GT . 100 ) GO  TO  120  SKI 

IGOBAK-O  SKI 

NELM=20 .00 /ELSTEP +1.00  SKI 

FR=  FRR  SKI 

MR=MRR  SKI 

00  10  IEL=1, NELH  SKI 

EL  = (I  EL  -1 ) * ELST  EP  SKI 

CALL  TRISL(R1, NL1, WL1,AZ,  EL , 0 . DO , FS , MODE, RS, RS 2( IEL > ,NLT »WLT  , SKI 

1PP  ,GP ( IEL)  ,AZF, ELF, SRANGE ,ELSR, D1ST, BEAR, IRETRN)  SKI 

IF(IRETRN-i)  6,8,12  SKI 

6  DO  7 1=1, INUSE  SKI 

IF(NNUSE(I ) .NE.3)  GO  TO  8 SKI 

7  CONTINUE  SKI 

GO  TO  10  SKI 

8  RS2  (I  EL  )=  1.020  SKI 

GP(IFL) =1.020  SKI 

10  CONTINUE  SKI 

IEL  = NEL  M+ 1 SKI 

12  IF(IEL.EQ.l)  GO  TO  110  SKI 

RS2 ( IEL ) = 1. 020  SKI 

NELM=IEL-1  SKI 

1 MI  N=  1 SKI 

IF=  0 SKI 

DO  20  IEL  = 1 * NELM  SKI 

I F( RS2 ( IEL)  ,GT. 6378. 8500)  GO  TO  20  SKI 

IF=  1 SKI 

IF(GP(IEL).LT.GP(IMIN))  IMIN  = IEL  SKI 

20  CONTINUE  SKI 

IFdF.EQ.  0)  GO  TO  110  SKI 

GPS=  GP ( IM IN)  SKI 

IFdOPT  .EO.O)  GO  TO  130  SKI 

IF(RS2( IMIN+1) ,LT. 6378. 8500)  GO  TO  130  SKI 

IF(RS2(  IMIN+D.GE. 0.9D20)  GO  TO  130  SKI 

IF(GP (IMIN  + 1) .GT.GP(IMIN)  ) GO  TO  130  SKI 

GPS=GP(  IMIN)  +(GP  (IMIN  + 1)  -GP  (I MIN)  )M  6378.  B5O0-RS2  (IMIN))/  SKI 

1 (RS2(IMIN+1) -RS2(IMIN))  SKI 

GO  TO  130  SKI 

110  IGOBA;<=  1 SKI 

GO  TO  130  SKI 

1 26  IGOBAK=  2 SKI 

130  OUT (1)  = FR  SKI 

OUT  (2)  =MR  SKI 

OUT  (3)  = GPS  SKI 

WRITE (6 '140)  OUT , IGOBAK  SKI 

140  FORMAT ( l SK  IP( ,2F1 0 . 6, F 10 . 2, 15)  SKI 

RETURN  SKI 

END  SKI 
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ROUTINE  PEAK 


74/74  DPT  = 1 


FT N 4.5*414 


04/27 


SUBROUTINE  PE AK < FRR , MRR, GPP , IGOBAK ) PEA 

IMPLICIT  DOUBLE  PRECISION (A-H, 0-7 > PEA 

REAL  0UT( 3 ) PEA 

DOUBLE  PRECISION  MRR, NL1 , NLT ,E LP f 2 5 ) , NLREr , MR  PEA 

COMMON/PERTC/NLREF,  W.  REF  , FR,  FO , MR,  DU  M<  1 1)  PEA 

COMMON/ FIOOLE/FP.FS ,P  1,NL 1 ,WL 1 , AZ , RS , EL  ST FP , MODE ,NS KPCL  PEA 

COMMON/ RE MWEN/INUSE,NNUSE( 20 ), FCREF ( 20>  PEA 

DATA  EL  P/ 11.900, 10.300, 8.900,,  7. 700,6. 7D0»5. 900, 5.2D0, 4.600,  PEA 

1 4. ID  0,3. 900, 5. 200, 4.  80 0,4. 500, 4. 200,4.00, 3. 800, 2*3.700,  PEA 

2 2*3.600,3*3,500,2*3.400/  PEA 

FR=FRR  PEA 

HR=  MRR  PEA 

IF=FP-4.74D0  PEA 

CALL  TRISL(R1,NL1,WL1,AZ,ELP(IF)  ,0.00, ^P, MODE, RS,  PEA 

1RS2,NLT,WLT,PP,GPP,AZF,ELF,SRANGE,ELSR,0IST, BEAR, IGOBAK)  PEA 

IF1RS2.GT. 6378. 8500)  IGOBAK=l  PEA 

DO  1 1=1, INUSE  PEA 

IF  (NNUSE(I ) .NE.3)  IG08AK  = 1 PEA 

1 CONTINUE  PEA 

OUT  <1)  = FR  PEA 

OUT 12)=MR  PEA 

OUT (3)  = GPP  PEA 

WRITE  (6,11)  OUT, IGOBAK  PEA 

11  FORMATU  PEAK!  ,2F10.6,F10.2,I5)  PEA 

RETURN  PEA 

END  PEA 


